mirror of
https://github.com/pfloos/quack
synced 2025-05-06 23:34:42 +02:00
complex diis and cap integrals
This commit is contained in:
Loris Burth
parent
e7c837ada7
commit
211168ecb9
213
PyDuck.py
Normal file → Executable file
213
PyDuck.py
Normal file → Executable file
@ -1,6 +1,7 @@
|
||||
#!/usr/bin/env python3
|
||||
|
||||
import os, sys
|
||||
import os
|
||||
import sys
|
||||
import argparse
|
||||
import pyscf
|
||||
from pyscf import gto
|
||||
@ -9,47 +10,60 @@ import subprocess
|
||||
import time
|
||||
|
||||
|
||||
#Find the value of the environnement variable QUACK_ROOT. If not present we use the current repository
|
||||
# Find the value of the environnement variable QUACK_ROOT. If not present we use the current repository
|
||||
if "QUACK_ROOT" not in os.environ:
|
||||
print("Please set the QUACK_ROOT environment variable, for example:\n")
|
||||
print("$ export QUACK_ROOT={0}".format(os.getcwd()))
|
||||
sys.exit(1)
|
||||
QuAcK_dir=os.environ.get('QUACK_ROOT','./')
|
||||
print("Please set the QUACK_ROOT environment variable, for example:\n")
|
||||
print("$ export QUACK_ROOT={0}".format(os.getcwd()))
|
||||
sys.exit(1)
|
||||
QuAcK_dir = os.environ.get('QUACK_ROOT', './')
|
||||
|
||||
#Create the argument parser object and gives a description of the script
|
||||
parser = argparse.ArgumentParser(description='This script is the main script of QuAcK, it is used to run the calculation.\n If $QUACK_ROOT is not set, $QUACK_ROOT is replaces by the current directory.')
|
||||
# Create the argument parser object and gives a description of the script
|
||||
parser = argparse.ArgumentParser(
|
||||
description='This script is the main script of QuAcK, it is used to run the calculation.\n If $QUACK_ROOT is not set, $QUACK_ROOT is replaces by the current directory.')
|
||||
|
||||
#Initialize all the options for the script
|
||||
parser.add_argument('-b', '--basis', type=str, required=True, help='Name of the file containing the basis set in the $QUACK_ROOT/basis/ directory')
|
||||
parser.add_argument('--bohr', default='Angstrom', action='store_const', const='Bohr', help='By default QuAcK assumes that the xyz files are in Angstrom. Add this argument if your xyz file is in Bohr.')
|
||||
parser.add_argument('-c', '--charge', type=int, default=0, help='Total charge of the molecule. Specify negative charges with "m" instead of the minus sign, for example m1 instead of -1. Default is 0')
|
||||
parser.add_argument('--cartesian', default=False, action='store_true', help='Add this option if you want to use cartesian basis functions.')
|
||||
parser.add_argument('--print_2e', default=True, action='store_true', help='If True, print 2e-integrals to disk.')
|
||||
parser.add_argument('--formatted_2e', default=False, action='store_true', help='Add this option if you want to print formatted 2e-integrals.')
|
||||
parser.add_argument('--mmap_2e', default=False, action='store_true', help='If True, avoid using DRAM when generating 2e-integrals.')
|
||||
parser.add_argument('--aosym_2e', default=False, action='store_true', help='If True, use 8-fold symmetry 2e-integrals.')
|
||||
parser.add_argument('-fc', '--frozen_core', type=bool, default=False, help='Freeze core MOs. Default is false')
|
||||
parser.add_argument('-m', '--multiplicity', type=int, default=1, help='Spin multiplicity. Default is 1 therefore singlet')
|
||||
parser.add_argument('--working_dir', type=str, default=QuAcK_dir, help='Set a working directory to run the calculation.')
|
||||
parser.add_argument('-x', '--xyz', type=str, required=True, help='Name of the file containing the nuclear coordinates in xyz format in the $QUACK_ROOT/mol/ directory without the .xyz extension')
|
||||
# Initialize all the options for the script
|
||||
parser.add_argument('-b', '--basis', type=str, required=True,
|
||||
help='Name of the file containing the basis set in the $QUACK_ROOT/basis/ directory')
|
||||
parser.add_argument('--bohr', default='Angstrom', action='store_const', const='Bohr',
|
||||
help='By default QuAcK assumes that the xyz files are in Angstrom. Add this argument if your xyz file is in Bohr.')
|
||||
parser.add_argument('-c', '--charge', type=int, default=0,
|
||||
help='Total charge of the molecule. Specify negative charges with "m" instead of the minus sign, for example m1 instead of -1. Default is 0')
|
||||
parser.add_argument('--cartesian', default=False, action='store_true',
|
||||
help='Add this option if you want to use cartesian basis functions.')
|
||||
parser.add_argument('--print_2e', default=True, action='store_true',
|
||||
help='If True, print 2e-integrals to disk.')
|
||||
parser.add_argument('--formatted_2e', default=False, action='store_true',
|
||||
help='Add this option if you want to print formatted 2e-integrals.')
|
||||
parser.add_argument('--mmap_2e', default=False, action='store_true',
|
||||
help='If True, avoid using DRAM when generating 2e-integrals.')
|
||||
parser.add_argument('--aosym_2e', default=False, action='store_true',
|
||||
help='If True, use 8-fold symmetry 2e-integrals.')
|
||||
parser.add_argument('-fc', '--frozen_core', type=bool,
|
||||
default=False, help='Freeze core MOs. Default is false')
|
||||
parser.add_argument('-m', '--multiplicity', type=int, default=1,
|
||||
help='Spin multiplicity. Default is 1 therefore singlet')
|
||||
parser.add_argument('--working_dir', type=str, default=QuAcK_dir,
|
||||
help='Set a working directory to run the calculation.')
|
||||
parser.add_argument('-x', '--xyz', type=str, required=True,
|
||||
help='Name of the file containing the nuclear coordinates in xyz format in the $QUACK_ROOT/mol/ directory without the .xyz extension')
|
||||
|
||||
#Parse the arguments
|
||||
# Parse the arguments
|
||||
args = parser.parse_args()
|
||||
input_basis=args.basis
|
||||
unit=args.bohr
|
||||
charge=args.charge
|
||||
frozen_core=args.frozen_core
|
||||
multiplicity=args.multiplicity
|
||||
xyz=args.xyz + '.xyz'
|
||||
cartesian=args.cartesian
|
||||
print_2e=args.print_2e
|
||||
formatted_2e=args.formatted_2e
|
||||
mmap_2e=args.mmap_2e
|
||||
aosym_2e=args.aosym_2e
|
||||
working_dir=args.working_dir
|
||||
input_basis = args.basis
|
||||
unit = args.bohr
|
||||
charge = args.charge
|
||||
frozen_core = args.frozen_core
|
||||
multiplicity = args.multiplicity
|
||||
xyz = args.xyz + '.xyz'
|
||||
cartesian = args.cartesian
|
||||
print_2e = args.print_2e
|
||||
formatted_2e = args.formatted_2e
|
||||
mmap_2e = args.mmap_2e
|
||||
aosym_2e = args.aosym_2e
|
||||
working_dir = args.working_dir
|
||||
|
||||
#Read molecule
|
||||
f = open(working_dir+'/mol/'+xyz,'r')
|
||||
# Read molecule
|
||||
f = open(working_dir+'/mol/'+xyz, 'r')
|
||||
lines = f.read().splitlines()
|
||||
nbAt = int(lines.pop(0))
|
||||
lines.pop(0)
|
||||
@ -57,98 +71,108 @@ list_pos_atom = []
|
||||
for line in lines:
|
||||
tmp = line.split()
|
||||
atom = tmp[0]
|
||||
pos = (float(tmp[1]),float(tmp[2]),float(tmp[3]))
|
||||
list_pos_atom.append([atom,pos])
|
||||
pos = (float(tmp[1]), float(tmp[2]), float(tmp[3]))
|
||||
list_pos_atom.append([atom, pos])
|
||||
f.close()
|
||||
|
||||
#Definition of the molecule
|
||||
# Definition of the molecule
|
||||
mol = gto.M(
|
||||
atom = list_pos_atom,
|
||||
basis = input_basis,
|
||||
charge = charge,
|
||||
spin = multiplicity - 1
|
||||
# symmetry = True # Enable symmetry
|
||||
atom=list_pos_atom,
|
||||
basis=input_basis,
|
||||
charge=charge,
|
||||
spin=multiplicity - 1
|
||||
# symmetry = True # Enable symmetry
|
||||
)
|
||||
|
||||
#Fix the unit for the lengths
|
||||
mol.unit=unit
|
||||
# Fix the unit for the lengths
|
||||
mol.unit = unit
|
||||
#
|
||||
mol.cart = cartesian
|
||||
|
||||
#Update mol object
|
||||
# Update mol object
|
||||
mol.build()
|
||||
|
||||
#Accessing number of electrons
|
||||
nelec=mol.nelec #Access the number of electrons
|
||||
nalpha=nelec[0]
|
||||
nbeta=nelec[1]
|
||||
# Accessing number of electrons
|
||||
nelec = mol.nelec # Access the number of electrons
|
||||
nalpha = nelec[0]
|
||||
nbeta = nelec[1]
|
||||
|
||||
subprocess.call(['mkdir', '-p', working_dir+'/input'])
|
||||
f = open(working_dir+'/input/molecule','w')
|
||||
f = open(working_dir+'/input/molecule', 'w')
|
||||
f.write('# nAt nEla nElb nCore nRyd\n')
|
||||
f.write(str(mol.natm)+' '+str(nalpha)+' '+str(nbeta)+' '+str(0)+' '+str(0)+'\n')
|
||||
f.write(str(mol.natm)+' '+str(nalpha)+' ' +
|
||||
str(nbeta)+' '+str(0)+' '+str(0)+'\n')
|
||||
f.write('# Znuc x y z\n')
|
||||
for i in range(len(list_pos_atom)):
|
||||
f.write(list_pos_atom[i][0]+' '+str(list_pos_atom[i][1][0])+' '+str(list_pos_atom[i][1][1])+' '+str(list_pos_atom[i][1][2])+'\n')
|
||||
f.write(list_pos_atom[i][0]+' '+str(list_pos_atom[i][1][0])+' ' +
|
||||
str(list_pos_atom[i][1][1])+' '+str(list_pos_atom[i][1][2])+'\n')
|
||||
f.close()
|
||||
|
||||
#Compute nuclear energy and put it in a file
|
||||
# Compute nuclear energy and put it in a file
|
||||
subprocess.call(['mkdir', '-p', working_dir+'/int'])
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/ENuc.dat'])
|
||||
f = open(working_dir+'/int/ENuc.dat','w')
|
||||
f = open(working_dir+'/int/ENuc.dat', 'w')
|
||||
f.write(str(mol.energy_nuc()))
|
||||
f.write(' ')
|
||||
f.close()
|
||||
|
||||
#Compute 1e integrals
|
||||
ovlp = mol.intor('int1e_ovlp')#Overlap matrix elements
|
||||
v1e = mol.intor('int1e_nuc') #Nuclear repulsion matrix elements
|
||||
t1e = mol.intor('int1e_kin') #Kinetic energy matrix elements
|
||||
dipole = mol.intor('int1e_r') #Matrix elements of the x, y, z operators
|
||||
x,y,z = dipole[0],dipole[1],dipole[2]
|
||||
# Compute 1e integrals
|
||||
ovlp = mol.intor('int1e_ovlp') # Overlap matrix elements
|
||||
v1e = mol.intor('int1e_nuc') # Nuclear repulsion matrix elements
|
||||
t1e = mol.intor('int1e_kin') # Kinetic energy matrix elements
|
||||
dipole = mol.intor('int1e_r') # Matrix elements of the x, y, z operators
|
||||
x, y, z = dipole[0], dipole[1], dipole[2]
|
||||
print(mol.nao)
|
||||
|
||||
norb = len(ovlp) # nBAS_AOs
|
||||
norb = len(ovlp) # nBAS_AOs
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/nBas.dat'])
|
||||
f = open(working_dir+'/int/nBas.dat','w')
|
||||
f = open(working_dir+'/int/nBas.dat', 'w')
|
||||
f.write(" {} ".format(str(norb)))
|
||||
f.close()
|
||||
|
||||
|
||||
def write_matrix_to_file(matrix,size,file,cutoff=1e-15):
|
||||
def write_matrix_to_file(matrix, size, file, cutoff=1e-15):
|
||||
f = open(file, 'w')
|
||||
for i in range(size):
|
||||
for j in range(i,size):
|
||||
for j in range(i, size):
|
||||
if abs(matrix[i][j]) > cutoff:
|
||||
f.write(str(i+1)+' '+str(j+1)+' '+"{:.16E}".format(matrix[i][j]))
|
||||
f.write(str(i+1)+' '+str(j+1)+' ' +
|
||||
"{:.16E}".format(matrix[i][j]))
|
||||
f.write('\n')
|
||||
f.close()
|
||||
|
||||
#Write all 1 electron quantities in files
|
||||
#Ov,Nuc,Kin,x,y,z
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/Ov.dat'])
|
||||
write_matrix_to_file(ovlp,norb,working_dir+'/int/Ov.dat')
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/Nuc.dat'])
|
||||
write_matrix_to_file(v1e,norb,working_dir+'/int/Nuc.dat')
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/Kin.dat'])
|
||||
write_matrix_to_file(t1e,norb,working_dir+'/int/Kin.dat')
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/x.dat'])
|
||||
write_matrix_to_file(x,norb,working_dir+'/int/x.dat')
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/y.dat'])
|
||||
write_matrix_to_file(y,norb,working_dir+'/int/y.dat')
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/z.dat'])
|
||||
write_matrix_to_file(z,norb,working_dir+'/int/z.dat')
|
||||
|
||||
def write_tensor_to_file(tensor,size,file_name,cutoff=1e-15):
|
||||
|
||||
# Write all 1 electron quantities in files
|
||||
# Ov,Nuc,Kin,x,y,z
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/Ov.dat'])
|
||||
write_matrix_to_file(ovlp, norb, working_dir+'/int/Ov.dat')
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/Nuc.dat'])
|
||||
write_matrix_to_file(v1e, norb, working_dir+'/int/Nuc.dat')
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/Kin.dat'])
|
||||
write_matrix_to_file(t1e, norb, working_dir+'/int/Kin.dat')
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/x.dat'])
|
||||
write_matrix_to_file(x, norb, working_dir+'/int/x.dat')
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/y.dat'])
|
||||
write_matrix_to_file(y, norb, working_dir+'/int/y.dat')
|
||||
subprocess.call(['rm', '-f', working_dir + '/int/z.dat'])
|
||||
write_matrix_to_file(z, norb, working_dir+'/int/z.dat')
|
||||
subprocess.call(['cp', '-r', working_dir +
|
||||
f'/int/cap/{args.xyz}/{input_basis}/CAP.dat', working_dir + f'/int/CAP.dat'])
|
||||
|
||||
|
||||
def write_tensor_to_file(tensor, size, file_name, cutoff=1e-15):
|
||||
f = open(file_name, 'w')
|
||||
for i in range(size):
|
||||
for j in range(i,size):
|
||||
for k in range(i,size):
|
||||
for l in range(j,size):
|
||||
for j in range(i, size):
|
||||
for k in range(i, size):
|
||||
for l in range(j, size):
|
||||
if abs(tensor[i][k][j][l]) > cutoff:
|
||||
f.write(str(i+1)+' '+str(j+1)+' '+str(k+1)+' '+str(l+1)+' '+"{:.16E}".format(tensor[i][k][j][l]))
|
||||
f.write(str(i+1)+' '+str(j+1)+' '+str(k+1)+' ' +
|
||||
str(l+1)+' '+"{:.16E}".format(tensor[i][k][j][l]))
|
||||
f.write('\n')
|
||||
f.close()
|
||||
|
||||
|
||||
if print_2e:
|
||||
# Write two-electron integrals to HD
|
||||
ti_2e = time.time()
|
||||
@ -158,7 +182,7 @@ if print_2e:
|
||||
subprocess.call(['rm', '-f', output_file_path])
|
||||
eri_ao = mol.intor('int2e')
|
||||
write_tensor_to_file(eri_ao, norb, output_file_path)
|
||||
|
||||
|
||||
if aosym_2e:
|
||||
output_file_path = working_dir + '/int/ERI_chem.bin'
|
||||
subprocess.call(['rm', '-f', output_file_path])
|
||||
@ -172,25 +196,24 @@ if print_2e:
|
||||
if(mmap_2e):
|
||||
# avoid using DRAM
|
||||
eri_shape = (norb, norb, norb, norb)
|
||||
eri_mmap = np.memmap(output_file_path, dtype='float64', mode='w+', shape=eri_shape)
|
||||
eri_mmap = np.memmap(
|
||||
output_file_path, dtype='float64', mode='w+', shape=eri_shape)
|
||||
mol.intor('int2e', out=eri_mmap)
|
||||
for i in range(norb):
|
||||
eri_mmap[i, :, :, :] = eri_mmap[i, :, :, :].transpose(1, 0, 2)
|
||||
eri_mmap.flush()
|
||||
del eri_mmap
|
||||
else:
|
||||
eri_ao = mol.intor('int2e').transpose(0, 2, 1, 3) # chem -> phys
|
||||
eri_ao = mol.intor('int2e').transpose(0, 2, 1, 3) # chem -> phys
|
||||
f = open(output_file_path, 'w')
|
||||
eri_ao.tofile(output_file_path)
|
||||
f.close()
|
||||
|
||||
te_2e = time.time()
|
||||
print("Wall time for writing 2e-integrals to disk: {:.3f} seconds".format(te_2e - ti_2e))
|
||||
print(
|
||||
"Wall time for writing 2e-integrals to disk: {:.3f} seconds".format(te_2e - ti_2e))
|
||||
sys.stdout.flush()
|
||||
|
||||
|
||||
|
||||
|
||||
#Execute the QuAcK fortran program
|
||||
# Execute the QuAcK fortran program
|
||||
subprocess.call([QuAcK_dir + '/bin/QuAcK', working_dir])
|
||||
|
||||
|
||||
@ -1,299 +0,0 @@
|
||||
#!/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")
|
||||
print(file_path)
|
||||
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_gram_schmidt(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 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 = orthonormalize(vecs)
|
||||
return vals, vecs
|
||||
|
||||
|
||||
def orthonormalize(vecs):
|
||||
# Orthonormalize them wrt cTc inner product
|
||||
R = vecs.T@vecs
|
||||
L = cholesky_decomposition(R)
|
||||
Linv = np.linalg.inv(L)
|
||||
vecs = vecs@Linv.T
|
||||
return 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)
|
||||
|
||||
|
||||
def DIIS_extrapolation(rcond, n_diis, error, e, error_in, e_inout):
|
||||
"""
|
||||
Perform DIIS extrapolation.
|
||||
|
||||
"""
|
||||
|
||||
# Update DIIS history by prepending new error and solution vectors
|
||||
error = np.column_stack((error_in, error[:, :-1])) # Shift history
|
||||
e = np.column_stack((e_inout, e[:, :-1])) # Shift history
|
||||
|
||||
# Build A matrix
|
||||
A = np.zeros((n_diis + 1, n_diis + 1), dtype=np.complex128)
|
||||
print(np.shape(error))
|
||||
A[:n_diis, :n_diis] = error@error.T
|
||||
A[:n_diis, n_diis] = -1.0
|
||||
A[n_diis, :n_diis] = -1.0
|
||||
A[n_diis, n_diis] = 0.0
|
||||
|
||||
# Build b vector
|
||||
b = np.zeros(n_diis + 1, dtype=np.complex128)
|
||||
b[n_diis] = -1.0
|
||||
|
||||
# Solve the linear system A * w = b
|
||||
try:
|
||||
w = np.linalg.solve(A, b)
|
||||
rcond = np.linalg.cond(A)
|
||||
except np.linalg.LinAlgError:
|
||||
raise ValueError("DIIS linear system is singular or ill-conditioned.")
|
||||
|
||||
# Extrapolate new solution
|
||||
e_inout[:] = w[:n_diis]@e[:, :n_diis].T
|
||||
|
||||
return rcond, n_diis, e_inout
|
||||
|
||||
|
||||
def cholesky_decomposition(A):
|
||||
"""
|
||||
Performs Cholesky-Decomposition wrt the c product. Returns L such that A = LTL
|
||||
"""
|
||||
|
||||
L = np.zeros_like(A)
|
||||
n = np.shape(L)[0]
|
||||
for i in range(n):
|
||||
for j in range(i + 1):
|
||||
s = A[i, j]
|
||||
|
||||
for k in range(j):
|
||||
s -= L[i, k] * L[j, k]
|
||||
|
||||
if i > j: # Off-diagonal elements
|
||||
L[i, j] = s / L[j, j]
|
||||
else: # Diagonal elements
|
||||
L[i, i] = s**0.5
|
||||
|
||||
return L
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
# Inputs
|
||||
workdir = "../"
|
||||
eta = 0.01
|
||||
thresh = 0.00001
|
||||
maxSCF = 256
|
||||
nO = read_nO(workdir)
|
||||
maxDiis = 5
|
||||
|
||||
# Read integrals
|
||||
T = read_matrix(f"{workdir}/int/Kin.dat")
|
||||
S = read_matrix(f"{workdir}/int/Ov.dat")
|
||||
V = read_matrix(f"{workdir}/int/Nuc.dat")
|
||||
ENuc = read_ENuc(f"{workdir}/int/ENuc.dat")
|
||||
nBas = np.shape(T)[0]
|
||||
nBasSq = nBas*nBas
|
||||
W = read_CAP_integrals(f"{workdir}/int/CAP.dat", nBas)
|
||||
ERI = read_2e_integrals(f"{workdir}/int/ERI.bin", nBas)
|
||||
X = get_X(S)
|
||||
W = -eta*W
|
||||
Hc = T + V + W*1j
|
||||
|
||||
# Initialization
|
||||
F_diis = np.zeros((nBasSq, maxDiis))
|
||||
error_diis = np.zeros((nBasSq, maxDiis))
|
||||
rcond = 0
|
||||
|
||||
# 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
|
||||
n_diis = 0
|
||||
|
||||
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
|
||||
|
||||
# # DIIS
|
||||
# n_diis = np.min([n_diis+1, maxDiis])
|
||||
# rcond, n_diis, F = DIIS_extrapolation(
|
||||
# rcond, n_diis, error_diis, F_diis, err, F)
|
||||
|
||||
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,41 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
|
||||
def generate_cap_file(filename, size, width, seed=42):
|
||||
"""
|
||||
Generates a file with random CAP integral values in the format:
|
||||
mu nu wx wy wz
|
||||
"""
|
||||
np.random.seed(seed) # For reproducibility
|
||||
with open(filename, 'w') as f:
|
||||
for mu in range(1, size + 1):
|
||||
for nu in range(mu, size + 1): # Only upper triangle to avoid duplicate entries
|
||||
# Generate three random values
|
||||
wx, wy, wz = np.random.rand(3)*width
|
||||
f.write(f"{mu} {nu} {wx:.6f} {wy:.6f} {wz:.6f}\n")
|
||||
|
||||
|
||||
def read_and_construct_matrix(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
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
with open("../int/nBas.dat", 'r') as f:
|
||||
size = int(f.readline().strip())
|
||||
print("nBas: ", size)
|
||||
width = 0
|
||||
generate_cap_file("../int/CAP.dat", size, width)
|
||||
W = read_and_construct_matrix("../int/CAP.dat", size)
|
||||
print("W matrix:")
|
||||
print(W)
|
||||
@ -2,5 +2,7 @@
|
||||
|
||||
cp ./methods.test ../input/methods
|
||||
cp ./options.test ../input/options
|
||||
basis=$2
|
||||
molecule=$1
|
||||
cd ..
|
||||
python3 PyDuck.py -x N2 -b sto-3g -c 0 -m 1
|
||||
python3 PyDuck.py -x $molecule -b $basis -c 0 -m 1
|
||||
|
||||
100
int/cap/CO/sto-3g/CAP.dat
Normal file
100
int/cap/CO/sto-3g/CAP.dat
Normal file
@ -0,0 +1,100 @@
|
||||
1 1 3.2591205169953606E-011 3.2591205169953606E-011 4.7270713420683057E-011
|
||||
1 2 1.0909850515273124E-007 1.0909850515273124E-007 1.0116692210295497E-007
|
||||
1 3 2.1201521662300532E-019 0.0000000000000000 0.0000000000000000
|
||||
1 4 0.0000000000000000 2.1201521662300532E-019 0.0000000000000000
|
||||
1 5 0.0000000000000000 0.0000000000000000 -1.7119541264283009E-007
|
||||
1 6 8.5110387749371516E-019 8.5110387749371525E-019 0.0000000000000000
|
||||
1 7 1.6053477154521912E-008 1.6053477154521912E-008 8.9580249908572546E-010
|
||||
1 8 1.5401421721512105E-020 0.0000000000000000 0.0000000000000000
|
||||
1 9 0.0000000000000000 1.5401421721512105E-020 0.0000000000000000
|
||||
1 10 -2.2009317144989998E-008 -2.2009317144989998E-008 -2.3185852966924291E-009
|
||||
2 1 1.0909850515273123E-007 1.0909850515273123E-007 1.0116692210295497E-007
|
||||
2 2 3.7657421170439367E-002 3.7657421170439367E-002 2.1213361581611412E-002
|
||||
2 3 8.6122699129564506E-018 0.0000000000000000 0.0000000000000000
|
||||
2 4 0.0000000000000000 8.6122699129564506E-018 0.0000000000000000
|
||||
2 5 0.0000000000000000 0.0000000000000000 -3.2960792036305107E-002
|
||||
2 6 1.1734955572834136E-011 1.1734955572834136E-011 2.7208889007654137E-015
|
||||
2 7 7.5133483187128083E-003 7.5133483187128074E-003 1.4697998957102060E-004
|
||||
2 8 6.4690907038125869E-018 0.0000000000000000 0.0000000000000000
|
||||
2 9 0.0000000000000000 6.4690907038125877E-018 0.0000000000000000
|
||||
2 10 -4.2626540158256449E-003 -4.2626540158256449E-003 6.2133296780354831E-005
|
||||
3 1 2.1201521662300535E-019 0.0000000000000000 0.0000000000000000
|
||||
3 2 8.6122699129564506E-018 0.0000000000000000 0.0000000000000000
|
||||
3 3 9.3785233366826731E-002 1.3493509930962135E-002 7.6845704009199586E-003
|
||||
3 4 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
3 5 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
3 6 3.3812080981053642E-021 0.0000000000000000 0.0000000000000000
|
||||
3 7 5.6846106957560671E-018 0.0000000000000000 0.0000000000000000
|
||||
3 8 1.8707831680260030E-002 2.4701269284935737E-003 4.6982188661143713E-005
|
||||
3 9 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
3 10 -4.2132123139296728E-018 0.0000000000000000 0.0000000000000000
|
||||
4 1 0.0000000000000000 2.1201521662300535E-019 0.0000000000000000
|
||||
4 2 0.0000000000000000 8.6122699129564506E-018 0.0000000000000000
|
||||
4 3 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
4 4 1.3493509930962135E-002 9.3785233366826731E-002 7.6845704009199586E-003
|
||||
4 5 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
4 6 0.0000000000000000 3.3812080981053646E-021 0.0000000000000000
|
||||
4 7 0.0000000000000000 5.6846106957560671E-018 0.0000000000000000
|
||||
4 8 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
4 9 2.4701269284935737E-003 1.8707831680260030E-002 4.6982188661143719E-005
|
||||
4 10 0.0000000000000000 -4.2132123139296728E-018 0.0000000000000000
|
||||
5 1 0.0000000000000000 0.0000000000000000 -1.7119541264283012E-007
|
||||
5 2 0.0000000000000000 0.0000000000000000 -3.2960792036305107E-002
|
||||
5 3 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
5 4 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
5 5 1.3493509930962135E-002 1.3493509930962135E-002 5.2329526092499988E-002
|
||||
5 6 1.5315965154753426E-011 1.5315965154753426E-011 5.7125594838600581E-015
|
||||
5 7 5.7504284972469356E-003 5.7504284972469356E-003 1.9245974391226696E-004
|
||||
5 8 5.9587410338227130E-018 0.0000000000000000 0.0000000000000000
|
||||
5 9 0.0000000000000000 5.9587410338227146E-018 0.0000000000000000
|
||||
5 10 -9.7283769522043325E-004 -9.7283769522043303E-004 6.2724149899989203E-004
|
||||
6 1 8.5110387749371516E-019 8.5110387749371525E-019 0.0000000000000000
|
||||
6 2 1.1734955572834136E-011 1.1734955572834136E-011 2.7208889007654137E-015
|
||||
6 3 3.3812080981053642E-021 0.0000000000000000 0.0000000000000000
|
||||
6 4 0.0000000000000000 3.3812080981053646E-021 0.0000000000000000
|
||||
6 5 1.5315965154753426E-011 1.5315965154753426E-011 5.7125594838600581E-015
|
||||
6 6 3.1228039277205289E-017 3.1228039277205289E-017 1.1243634227890234E-022
|
||||
6 7 3.3271555055667935E-011 3.3271555055667935E-011 3.9492252102805315E-014
|
||||
6 8 4.9206623360812111E-020 0.0000000000000000 0.0000000000000000
|
||||
6 9 0.0000000000000000 4.9206623360812117E-020 0.0000000000000000
|
||||
6 10 0.0000000000000000 0.0000000000000000 6.4605147029897322E-014
|
||||
7 1 1.6053477154521912E-008 1.6053477154521912E-008 8.9580249908572546E-010
|
||||
7 2 7.5133483187128083E-003 7.5133483187128074E-003 1.4697998957102060E-004
|
||||
7 3 5.6846106957560671E-018 0.0000000000000000 0.0000000000000000
|
||||
7 4 0.0000000000000000 5.6846106957560671E-018 0.0000000000000000
|
||||
7 5 5.7504284972469356E-003 5.7504284972469356E-003 1.9245974391226696E-004
|
||||
7 6 3.3271555055667935E-011 3.3271555055667935E-011 3.9492252102805315E-014
|
||||
7 7 6.8439309685517672E-003 6.8439309685517672E-003 1.3388069661601368E-003
|
||||
7 8 6.5258921023947781E-018 0.0000000000000000 0.0000000000000000
|
||||
7 9 0.0000000000000000 6.5258921023947781E-018 0.0000000000000000
|
||||
7 10 2.4002938818070129E-019 2.4002938818070129E-019 2.3454583731694427E-003
|
||||
8 1 1.5401421721512105E-020 0.0000000000000000 0.0000000000000000
|
||||
8 2 6.4690907038125869E-018 0.0000000000000000 0.0000000000000000
|
||||
8 3 1.8707831680260027E-002 2.4701269284935737E-003 4.6982188661143713E-005
|
||||
8 4 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
8 5 5.9587410338227130E-018 0.0000000000000000 0.0000000000000000
|
||||
8 6 4.9206623360812111E-020 0.0000000000000000 0.0000000000000000
|
||||
8 7 6.5258921023947781E-018 0.0000000000000000 0.0000000000000000
|
||||
8 8 1.9123613468117230E-002 2.2741125775108282E-003 4.3174144495353406E-004
|
||||
8 9 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
8 10 1.7559510230176076E-033 0.0000000000000000 0.0000000000000000
|
||||
9 1 0.0000000000000000 1.5401421721512105E-020 0.0000000000000000
|
||||
9 2 0.0000000000000000 6.4690907038125877E-018 0.0000000000000000
|
||||
9 3 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
9 4 2.4701269284935737E-003 1.8707831680260030E-002 4.6982188661143719E-005
|
||||
9 5 0.0000000000000000 5.9587410338227146E-018 0.0000000000000000
|
||||
9 6 0.0000000000000000 4.9206623360812117E-020 0.0000000000000000
|
||||
9 7 0.0000000000000000 6.5258921023947781E-018 0.0000000000000000
|
||||
9 8 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
9 9 2.2741125775108291E-003 1.9123613468117237E-002 4.3174144495353422E-004
|
||||
9 10 0.0000000000000000 1.7559510230176079E-033 0.0000000000000000
|
||||
10 1 -2.2009317144989998E-008 -2.2009317144989998E-008 -2.3185852966924291E-009
|
||||
10 2 -4.2626540158256449E-003 -4.2626540158256449E-003 6.2133296780354845E-005
|
||||
10 3 -4.2132123139296720E-018 0.0000000000000000 0.0000000000000000
|
||||
10 4 0.0000000000000000 -4.2132123139296720E-018 0.0000000000000000
|
||||
10 5 -9.7283769522043336E-004 -9.7283769522043303E-004 6.2724149899989203E-004
|
||||
10 6 0.0000000000000000 0.0000000000000000 6.4605147029897322E-014
|
||||
10 7 2.4002938818070129E-019 2.4002938818070129E-019 2.3454583731694427E-003
|
||||
10 8 1.7559510230176076E-033 0.0000000000000000 0.0000000000000000
|
||||
10 9 0.0000000000000000 1.7559510230176079E-033 0.0000000000000000
|
||||
10 10 2.2741125775108287E-003 2.2741125775108291E-003 4.1770865567530291E-003
|
||||
12100
int/cap/N2/aug-cc-pvtz/CAP.dat
Normal file
12100
int/cap/N2/aug-cc-pvtz/CAP.dat
Normal file
File diff suppressed because it is too large
Load Diff
4900
int/cap/N2/cc-pvtz/CAP.dat
Normal file
4900
int/cap/N2/cc-pvtz/CAP.dat
Normal file
File diff suppressed because it is too large
Load Diff
100
int/cap/N2/sto-3g/CAP.dat
Normal file
100
int/cap/N2/sto-3g/CAP.dat
Normal file
@ -0,0 +1,100 @@
|
||||
1 1 4.6655185822660997E-014 4.6655185822660997E-014 1.6497753194259012E-015
|
||||
1 2 2.2954533580438234E-009 2.2954533580438234E-009 2.5724228932520038E-010
|
||||
1 3 1.0127179643257497E-019 0.0000000000000000 0.0000000000000000
|
||||
1 4 0.0000000000000000 1.0127179643257497E-019 0.0000000000000000
|
||||
1 5 0.0000000000000000 0.0000000000000000 -4.3504742730995192E-010
|
||||
1 6 1.2747662113999374E-018 1.2747662113999376E-018 0.0000000000000000
|
||||
1 7 5.9476138926704675E-010 5.9476138926704664E-010 8.2273863063969177E-012
|
||||
1 8 8.0628842609266262E-021 0.0000000000000000 0.0000000000000000
|
||||
1 9 0.0000000000000000 8.0628842609266262E-021 0.0000000000000000
|
||||
1 10 -7.6854936906706609E-010 -7.6854936906706609E-010 -1.8607731397469347E-011
|
||||
2 1 2.2954533580438229E-009 2.2954533580438229E-009 2.5724228932520038E-010
|
||||
2 2 1.7873802543636113E-002 1.7873802543636113E-002 6.9594592588476199E-003
|
||||
2 3 1.5868624574356713E-017 0.0000000000000000 0.0000000000000000
|
||||
2 4 0.0000000000000000 1.5868624574356713E-017 0.0000000000000000
|
||||
2 5 0.0000000000000000 0.0000000000000000 -1.1268106770080917E-002
|
||||
2 6 5.9476138926704675E-010 5.9476138926704664E-010 8.2273863088512740E-012
|
||||
2 7 9.4886856622103648E-003 9.4886856622103648E-003 2.3568662049465330E-004
|
||||
2 8 7.0037535737045501E-018 0.0000000000000000 0.0000000000000000
|
||||
2 9 0.0000000000000000 7.0037535737045501E-018 0.0000000000000000
|
||||
2 10 -6.1569807285705019E-003 -6.1569807285705019E-003 -1.2521548512534853E-004
|
||||
3 1 1.0127179643257497E-019 0.0000000000000000 0.0000000000000000
|
||||
3 2 1.5868624574356713E-017 0.0000000000000000 0.0000000000000000
|
||||
3 3 4.6569688658962910E-002 6.1623613278307225E-003 2.3621922377360617E-003
|
||||
3 4 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
3 5 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
3 6 8.0628842609266277E-021 0.0000000000000000 0.0000000000000000
|
||||
3 7 7.0037535737045501E-018 0.0000000000000000 0.0000000000000000
|
||||
3 8 2.3466435264119647E-002 3.1721254678129711E-003 7.7747264619501400E-005
|
||||
3 9 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
3 10 -5.9826673506180821E-018 0.0000000000000000 0.0000000000000000
|
||||
4 1 0.0000000000000000 1.0127179643257497E-019 0.0000000000000000
|
||||
4 2 0.0000000000000000 1.5868624574356713E-017 0.0000000000000000
|
||||
4 3 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
4 4 6.1623613278307217E-003 4.6569688658962910E-002 2.3621922377360617E-003
|
||||
4 5 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
4 6 0.0000000000000000 8.0628842609266262E-021 0.0000000000000000
|
||||
4 7 0.0000000000000000 7.0037535737045501E-018 0.0000000000000000
|
||||
4 8 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
4 9 3.1721254678129707E-003 2.3466435264119644E-002 7.7747264619501400E-005
|
||||
4 10 0.0000000000000000 -5.9826673506180821E-018 0.0000000000000000
|
||||
5 1 0.0000000000000000 0.0000000000000000 -4.3504742730995192E-010
|
||||
5 2 0.0000000000000000 0.0000000000000000 -1.1268106770080917E-002
|
||||
5 3 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
5 4 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
5 5 6.1623613278307225E-003 6.1623613278307217E-003 1.8598627518531509E-002
|
||||
5 6 7.6854936906706609E-010 7.6854936906706609E-010 1.8607731405976441E-011
|
||||
5 7 6.1569807285705019E-003 6.1569807285705019E-003 1.2521548512535693E-004
|
||||
5 8 5.9826673506180821E-018 0.0000000000000000 0.0000000000000000
|
||||
5 9 0.0000000000000000 5.9826673506180821E-018 0.0000000000000000
|
||||
5 10 -9.6996539344591163E-004 -9.6996539344591163E-004 9.5475777343925674E-004
|
||||
6 1 1.2747662113999374E-018 1.2747662113999376E-018 0.0000000000000000
|
||||
6 2 5.9476138926704675E-010 5.9476138926704664E-010 8.2273863088512740E-012
|
||||
6 3 8.0628842609266262E-021 0.0000000000000000 0.0000000000000000
|
||||
6 4 0.0000000000000000 8.0628842609266262E-021 0.0000000000000000
|
||||
6 5 7.6854936906706609E-010 7.6854936906706609E-010 1.8607731405976441E-011
|
||||
6 6 4.6655185822660997E-014 4.6655185822660997E-014 1.6497780571212942E-015
|
||||
6 7 2.2954533580438234E-009 2.2954533580438234E-009 2.5724228934092053E-010
|
||||
6 8 1.0127179643257497E-019 0.0000000000000000 0.0000000000000000
|
||||
6 9 0.0000000000000000 1.0127179643257497E-019 0.0000000000000000
|
||||
6 10 0.0000000000000000 0.0000000000000000 4.3504742739585608E-010
|
||||
7 1 5.9476138926704675E-010 5.9476138926704664E-010 8.2273863063969177E-012
|
||||
7 2 9.4886856622103648E-003 9.4886856622103648E-003 2.3568662049465330E-004
|
||||
7 3 7.0037535737045501E-018 0.0000000000000000 0.0000000000000000
|
||||
7 4 0.0000000000000000 7.0037535737045501E-018 0.0000000000000000
|
||||
7 5 6.1569807285705019E-003 6.1569807285705019E-003 1.2521548512535696E-004
|
||||
7 6 2.2954533580438229E-009 2.2954533580438229E-009 2.5724228934092053E-010
|
||||
7 7 1.7873802543636113E-002 1.7873802543636113E-002 6.9594592588476216E-003
|
||||
7 8 1.5868624574356713E-017 0.0000000000000000 0.0000000000000000
|
||||
7 9 0.0000000000000000 1.5868624574356713E-017 0.0000000000000000
|
||||
7 10 -8.4737143932112045E-019 -8.4737143932112045E-019 1.1268106770080933E-002
|
||||
8 1 8.0628842609266277E-021 0.0000000000000000 0.0000000000000000
|
||||
8 2 7.0037535737045501E-018 0.0000000000000000 0.0000000000000000
|
||||
8 3 2.3466435264119647E-002 3.1721254678129711E-003 7.7747264619501400E-005
|
||||
8 4 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
8 5 5.9826673506180821E-018 0.0000000000000000 0.0000000000000000
|
||||
8 6 1.0127179643257497E-019 0.0000000000000000 0.0000000000000000
|
||||
8 7 1.5868624574356713E-017 0.0000000000000000 0.0000000000000000
|
||||
8 8 4.6569688658962910E-002 6.1623613278307225E-003 2.3621922377360612E-003
|
||||
8 9 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
8 10 -3.0645883512780251E-033 0.0000000000000000 0.0000000000000000
|
||||
9 1 0.0000000000000000 8.0628842609266262E-021 0.0000000000000000
|
||||
9 2 0.0000000000000000 7.0037535737045501E-018 0.0000000000000000
|
||||
9 3 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
9 4 3.1721254678129707E-003 2.3466435264119644E-002 7.7747264619501400E-005
|
||||
9 5 0.0000000000000000 5.9826673506180821E-018 0.0000000000000000
|
||||
9 6 0.0000000000000000 1.0127179643257497E-019 0.0000000000000000
|
||||
9 7 0.0000000000000000 1.5868624574356713E-017 0.0000000000000000
|
||||
9 8 0.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
9 9 6.1623613278307217E-003 4.6569688658962910E-002 2.3621922377360612E-003
|
||||
9 10 0.0000000000000000 -3.0645883512780251E-033 0.0000000000000000
|
||||
10 1 -7.6854936906706609E-010 -7.6854936906706609E-010 -1.8607731397469347E-011
|
||||
10 2 -6.1569807285705019E-003 -6.1569807285705019E-003 -1.2521548512534853E-004
|
||||
10 3 -5.9826673506180829E-018 0.0000000000000000 0.0000000000000000
|
||||
10 4 0.0000000000000000 -5.9826673506180829E-018 0.0000000000000000
|
||||
10 5 -9.6996539344591163E-004 -9.6996539344591163E-004 9.5475777343925685E-004
|
||||
10 6 0.0000000000000000 0.0000000000000000 4.3504742739585608E-010
|
||||
10 7 -8.4737143932112045E-019 -8.4737143932112045E-019 1.1268106770080933E-002
|
||||
10 8 -3.0645883512780251E-033 0.0000000000000000 0.0000000000000000
|
||||
10 9 0.0000000000000000 -3.0645883512780251E-033 0.0000000000000000
|
||||
10 10 6.1623613278307225E-003 6.1623613278307217E-003 1.8598627518531523E-002
|
||||
@ -1,5 +1,5 @@
|
||||
subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nO,S,T,V,ERI,dipole_int,X,ERHF,eHF,c,P)
|
||||
subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
|
||||
nBas,nO,S,T,V,ERI,X,ERHF,eHF,c,P,F)
|
||||
|
||||
! Perform complex restricted Hartree-Fock calculation
|
||||
|
||||
@ -18,17 +18,12 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nNuc
|
||||
double precision,intent(in) :: ZNuc(nNuc)
|
||||
double precision,intent(in) :: rNuc(nNuc,ncart)
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: S(nBas,nBas)
|
||||
double precision,intent(in) :: T(nBas,nBas)
|
||||
double precision,intent(in) :: V(nBas,nBas)
|
||||
double precision,intent(in) :: X(nBas,nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nSCF
|
||||
@ -38,7 +33,6 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
|
||||
complex*16 :: EV
|
||||
complex*16 :: EJ
|
||||
complex*16 :: EK
|
||||
complex*16 :: dipole(ncart)
|
||||
|
||||
double precision :: Conv
|
||||
double precision :: rcond
|
||||
@ -46,16 +40,15 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
|
||||
|
||||
double precision :: eta
|
||||
double precision,allocatable :: W(:,:)
|
||||
complex*16,allocatable :: Hc(:,:)
|
||||
complex*16,allocatable :: J(:,:)
|
||||
complex*16,allocatable :: K(:,:)
|
||||
complex*16,allocatable :: cp(:,:)
|
||||
complex*16,allocatable :: F(:,:)
|
||||
complex*16,allocatable :: Fp(:,:)
|
||||
complex*16,allocatable :: err(:,:)
|
||||
complex*16,allocatable :: err_diis(:,:)
|
||||
complex*16,allocatable :: F_diis(:,:)
|
||||
complex*16,allocatable :: Z(:,:)
|
||||
complex*16,allocatable :: Hc(:,:)
|
||||
|
||||
|
||||
! Output variables
|
||||
|
||||
@ -63,6 +56,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
|
||||
complex*16,intent(out) :: eHF(nBas)
|
||||
complex*16,intent(inout) :: c(nBas,nBas)
|
||||
complex*16,intent(out) :: P(nBas,nBas)
|
||||
complex*16,intent(inout) :: F(nBas,nBas)
|
||||
|
||||
! Hello world
|
||||
|
||||
@ -79,9 +73,16 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
|
||||
|
||||
! Memory allocation
|
||||
|
||||
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))
|
||||
allocate(err_diis(nBasSq,max_diis))
|
||||
allocate(F_diis(nBasSq,max_diis))
|
||||
allocate(Hc(nBas,nBas))
|
||||
allocate(W(nBas,nBas))
|
||||
allocate(J(nBas,nBas))
|
||||
allocate(K(nBas,nBas))
|
||||
allocate(err(nBas,nBas))
|
||||
allocate(cp(nBas,nBas))
|
||||
allocate(Fp(nBas,nBas))
|
||||
|
||||
! Read CAP integrals from file
|
||||
call read_CAP_integrals(nBas,W)
|
||||
W(:,:) = -eta*W(:,:)
|
||||
@ -97,25 +98,23 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
|
||||
n_diis = 0
|
||||
F_diis(:,:) = cmplx(0d0,0d0,kind=8)
|
||||
err_diis(:,:) = cmplx(0d0,0d0,kind=8)
|
||||
rcond = 0d0
|
||||
|
||||
rcond = 0d0
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'--------------------------------------------------------------------------------------------------'
|
||||
write(*,*)'-------------------------------------------------------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A36,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') &
|
||||
'|','#','|','E(RHF)','|','EJ(RHF)','|','EK(RHF)','|','Conv','|'
|
||||
write(*,*)'--------------------------------------------------------------------------------------------------'
|
||||
write(*,*)'-------------------------------------------------------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Build Fock matrix
|
||||
|
||||
call complex_Hartree_matrix_AO_basis(nBas,P,ERI,J)
|
||||
@ -150,18 +149,17 @@ 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 !
|
||||
|
||||
! 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)
|
||||
!
|
||||
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)
|
||||
if (rcond<1.0d-10) write(*,*) "!!! DIIS system ill conditioned: rcond = ", rcond," !!!"
|
||||
end if
|
||||
! 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(:,:)))
|
||||
@ -175,9 +173,9 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
|
||||
! Dump results
|
||||
|
||||
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,'|'
|
||||
'|',nSCF,'|',real(ERHF + ENuc),'+',aimag(ERHF),'i','|',real(EJ),'|',real(EK),'|',Conv,'|'
|
||||
end do
|
||||
write(*,*)'--------------------------------------------------------------------------------------------------'
|
||||
write(*,*)'-------------------------------------------------------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
@ -196,9 +194,9 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
|
||||
|
||||
end if
|
||||
|
||||
! Compute dipole moments
|
||||
call print_cRHF(nBas,nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF)
|
||||
! Testing zone
|
||||
|
||||
! Testing zone
|
||||
|
||||
if(dotest) then
|
||||
|
||||
@ -208,5 +206,5 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
|
||||
! call dump_test_value('R','RHF dipole moment',norm2(dipole))
|
||||
|
||||
end if
|
||||
|
||||
deallocate(J,K,err,cp,Fp,err_diis,F_diis,Hc)
|
||||
end subroutine
|
||||
|
||||
@ -23,7 +23,6 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: ixyz
|
||||
integer :: HOMO
|
||||
integer :: LUMO
|
||||
complex*16 :: Gap
|
||||
@ -43,40 +42,41 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF)
|
||||
! Dump results
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A50)') '------------------------------------------------------------'
|
||||
write(*,'(A69)') '---------------------------------------------------------'
|
||||
write(*,'(A33)') ' Summary '
|
||||
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,A1,F16.10,A1,A3)') ' cRHF energy = ',real(ERHF + ENuc),'+',aimag(ERHF+ENuc),'i',' au'
|
||||
write(*,'(A50)') '------------------------------------------------------------'
|
||||
write(*,'(A69)') '---------------------------------------------------------'
|
||||
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' One-electron energy = ',real(ET + EV),'+',aimag(ET+EV),'i',' au'
|
||||
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Kinetic energy = ',real(ET),'+',aimag(ET),'i',' au'
|
||||
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Potential energy = ',real(EV),'+',aimag(Ev),'i',' au'
|
||||
write(*,'(A69)') '---------------------------------------------------'
|
||||
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Two-electron energy = ',real(EJ + EK),'+',aimag(EJ+EK),'i',' au'
|
||||
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Hartree energy = ',real(EJ),'+',aimag(EJ),'i',' au'
|
||||
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Exchange energy = ',real(EK),'+',aimag(EK),'i',' au'
|
||||
write(*,'(A69)') '---------------------------------------------------------'
|
||||
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Electronic energy = ',real(ERHF),'+',aimag(ERHF),'i',' au'
|
||||
write(*,'(A33,1X,F16.10,19X,A3)') ' Nuclear repulsion = ',ENuc,' au'
|
||||
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' cRHF energy = ',real(ERHF + ENuc),'+',aimag(ERHF+ENuc),'i',' au'
|
||||
write(*,'(A69)') '---------------------------------------------------------'
|
||||
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(*,'(A69)') '---------------------------------------------------------'
|
||||
write(*,'(A33,1X,F16.6)') ' <Sz> = ',S
|
||||
write(*,'(A33,1X,F16.6)') ' <S^2> = ',S2
|
||||
write(*,*)
|
||||
|
||||
! Print results
|
||||
|
||||
if(dump_orb) then
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A50)') ' cRHF orbital coefficients '
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A69)') '---------------------------------------'
|
||||
write(*,'(A69)') ' cRHF orbital coefficients '
|
||||
write(*,'(A69)') '---------------------------------------'
|
||||
call complex_matout(nBas, nOrb, cHF)
|
||||
write(*,*)
|
||||
end if
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A50)') ' cRHF orbital energies (au) '
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A37)') '-------------------------------'
|
||||
write(*,'(A37)') ' cRHF orbital energies (au) '
|
||||
write(*,'(A37)') '-------------------------------'
|
||||
call complex_vecout(nOrb, eHF)
|
||||
write(*,*)
|
||||
|
||||
|
||||
@ -321,8 +321,11 @@ program QuAcK
|
||||
write(*,*)
|
||||
|
||||
! Memory deallocation
|
||||
|
||||
deallocate(ZNuc,rNuc)
|
||||
deallocate(S,T,V,Hc,dipole_int_AO)
|
||||
|
||||
if (allocated(rNuc)) deallocate(rNuc)
|
||||
if (allocated(Znuc)) deallocate(Znuc)
|
||||
if (allocated(T)) deallocate(T)
|
||||
if (allocated(V)) deallocate(V)
|
||||
if (allocated(Hc)) deallocate(Hc)
|
||||
if (allocated(dipole_int_AO)) deallocate(dipole_int_AO)
|
||||
!if (allocated(S)) deallocate(S)
|
||||
end program
|
||||
|
||||
@ -104,6 +104,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
|
||||
double precision,allocatable :: PHF(:,:)
|
||||
complex*16,allocatable :: complex_PHF(:,:)
|
||||
double precision,allocatable :: FHF(:,:)
|
||||
complex*16,allocatable :: complex_FHF(:,:)
|
||||
double precision :: ERHF
|
||||
complex*16 :: complex_ERHF
|
||||
double precision,allocatable :: dipole_int_MO(:,:,:)
|
||||
@ -127,6 +128,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
|
||||
allocate(complex_eHF(nOrb))
|
||||
allocate(complex_cHF(nBas,nOrb))
|
||||
allocate(complex_PHF(nBas,nBas))
|
||||
allocate(complex_FHF(nBas,nBas))
|
||||
allocate(ERI_AO(nBas,nBas,nBas,nBas))
|
||||
allocate(FHF(nBas,nBas))
|
||||
allocate(dipole_int_MO(nOrb,nOrb,ncart))
|
||||
@ -171,8 +173,8 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
|
||||
|
||||
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,complex_eHF,complex_cHF,complex_PHF)
|
||||
call cRHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,ENuc, &
|
||||
nBas,nO,S,T,V,ERI_AO,X,complex_ERHF,complex_eHF,complex_cHF,complex_PHF,complex_FHF)
|
||||
call wall_time(end_HF)
|
||||
|
||||
t_HF = end_HF - start_HF
|
||||
@ -390,4 +392,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
|
||||
deallocate(ERI_MO)
|
||||
deallocate(ERI_AO)
|
||||
|
||||
deallocate(complex_eHF)
|
||||
deallocate(complex_cHF)
|
||||
deallocate(complex_PHF)
|
||||
end subroutine
|
||||
|
||||
@ -29,10 +29,9 @@ subroutine complex_DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_
|
||||
allocate(A(n_diis+1,n_diis+1),b(n_diis+1),w(n_diis+1))
|
||||
|
||||
! Update DIIS "history"
|
||||
|
||||
call complex_prepend(n_err,n_diis,error,error_in)
|
||||
call complex_prepend(n_e,n_diis,e,e_inout)
|
||||
|
||||
|
||||
! Build A matrix
|
||||
|
||||
A(1:n_diis,1:n_diis) = matmul(transpose(error),error)
|
||||
@ -47,12 +46,10 @@ subroutine complex_DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_
|
||||
b(n_diis+1) = cmplx(-1d0,0d0,kind=8)
|
||||
|
||||
! Solve linear system
|
||||
call complex_vecout(b)
|
||||
call complex_matout(A)
|
||||
call complex_linear_solve(n_diis+1,A,b,w,rcond)
|
||||
|
||||
! Extrapolate
|
||||
|
||||
e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis)))
|
||||
|
||||
deallocate(A,b,w)
|
||||
end subroutine
|
||||
|
||||
@ -7,9 +7,9 @@ subroutine complex_matout(n,m,A)
|
||||
integer,intent(in) :: n,m
|
||||
complex*16,intent(in) :: A(n,m)
|
||||
|
||||
write( *,*) 'Real part:'
|
||||
write( *,*) 'Re'
|
||||
call matout(n,m,real(A))
|
||||
write (*,*) 'Imaginary part:'
|
||||
write (*,*) 'Im'
|
||||
call matout(n,m,aimag(A))
|
||||
|
||||
end subroutine
|
||||
|
||||
@ -9,7 +9,7 @@ subroutine complex_vecout(n,v)
|
||||
double precision,allocatable :: v2(:,:)
|
||||
|
||||
allocate(v2(n,2))
|
||||
write(*,*) 'First column real part, second imaginary part'
|
||||
write(*,'(17X,A2,13X,A2)') 'Re','Im'
|
||||
v2(:,1) = real(v)
|
||||
v2(:,2) = aimag(v)
|
||||
call matout(n,2,v2)
|
||||
|
||||
@ -305,7 +305,7 @@ subroutine complex_prepend(N,M,A,b)
|
||||
integer,intent(in) :: N,M
|
||||
complex*16,intent(in) :: b(N)
|
||||
|
||||
! Local viaruabkes
|
||||
! Local variables
|
||||
|
||||
integer :: i,j
|
||||
|
||||
@ -315,7 +315,7 @@ subroutine complex_prepend(N,M,A,b)
|
||||
|
||||
|
||||
! print*,'b in append'
|
||||
! call matout(N,1,b)
|
||||
! call complex_matout(N,1,b)
|
||||
|
||||
do i=1,N
|
||||
do j=M-1,1,-1
|
||||
|
||||
@ -334,31 +334,50 @@ subroutine linear_solve(N,A,b,x,rcond)
|
||||
! stop 'error in linear_solve (dsysvx)!!'
|
||||
|
||||
! end if
|
||||
|
||||
deallocate(work,ipiv,iwork,AF)
|
||||
end subroutine
|
||||
|
||||
subroutine complex_linear_solve(N,A,b,x,rcond)
|
||||
|
||||
! Solve the linear system A.x = b where A is a NxN matrix
|
||||
! and x and b are vectors of size N
|
||||
! and x and x are vectors of size N
|
||||
|
||||
implicit none
|
||||
|
||||
integer,intent(in) :: N
|
||||
complex*16,intent(out) :: A(N,N),b(N)
|
||||
double precision :: rcond
|
||||
complex*16,intent(out) :: x(N)
|
||||
double precision,intent(out) :: rcond
|
||||
|
||||
integer :: info
|
||||
integer :: info,lwork
|
||||
double precision :: ferr,berr
|
||||
integer,allocatable :: ipiv(:)
|
||||
double precision,allocatable :: rwork(:)
|
||||
complex*16,allocatable :: AF(:,:),work(:)
|
||||
|
||||
! Find optimal size for temporary arrays
|
||||
|
||||
allocate(ipiv(N))
|
||||
call zgesv(N,1,A,N,ipiv,b,N,info)
|
||||
allocate(work(1))
|
||||
allocate(AF(N,N),ipiv(N),rwork(N))
|
||||
|
||||
end subroutine
|
||||
lwork = -1
|
||||
call zsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,rwork,info)
|
||||
lwork = max(1,int(real(work(1))))
|
||||
|
||||
deallocate(work)
|
||||
|
||||
allocate(work(lwork))
|
||||
|
||||
call zsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,rwork,info)
|
||||
|
||||
if (info /= 0) then
|
||||
|
||||
print *, info
|
||||
stop 'error in linear_solve (zsysv)!!'
|
||||
|
||||
end if
|
||||
deallocate(work,ipiv,rwork,AF)
|
||||
end subroutine
|
||||
|
||||
subroutine easy_linear_solve(N,A,b,x)
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user