LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-05-06 23:34:42 +02:00
Code Releases Activity

G0W0 corrected with CAP

Browse Source
This commit is contained in:
Loris Burth 2025-03-18 15:46:32 +01:00
parent 8ab381ed27
commit d7418ae352
17 changed files with 205 additions and 523 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
PyDuck.py
View File

@ -8,6 +8,7 @@ from pyscf import gto
import numpy as np
import subprocess
import time
import pyopencap
# Find the value of the environnement variable QUACK_ROOT. If not present we use the current repository
@ -23,7 +24,9 @@ parser = argparse.ArgumentParser(
# 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')
help='Name of the file containing the basis set in the $QUACK_ROOT/basis/ directory. If cap is used the basis in psi4 style has to be provided in the directory cap_data/basis directory.')
parser.add_argument('--use_local_basis', default=False, action='store_true',
help='If True, basis is loaded from local storage. Needed for CAP.')
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,
@ -46,10 +49,15 @@ 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')
parser.add_argument("--use_cap", action="store_true", default=False,
help="If true cap integrals are calculated by opencap and written to a file. The basis has to be provided in the cap_data/basis dir in psi4 style and the onsets in the cap_data/onsets dir.")
# Parse the arguments
args = parser.parse_args()
input_basis = args.basis
# Basisset path
pyscf_path = os.path.dirname(pyscf.__file__)
basis_path = os.path.join(pyscf_path, "gto", "basis", input_basis + ".dat")
unit = args.bohr
charge = args.charge
frozen_core = args.frozen_core
@ -61,7 +69,6 @@ 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')
lines = f.read().splitlines()
@ -122,7 +129,6 @@ 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
subprocess.call(['rm', '-f', working_dir + '/int/nBas.dat'])
@ -130,6 +136,44 @@ f = open(working_dir+'/int/nBas.dat', 'w')
f.write(" {} ".format(str(norb)))
f.close()
# CAP definition
if args.use_cap:
f = open(working_dir+'/cap_data/onsets/'+args.xyz, 'r')
lines = f.read().splitlines()
for line in lines:
tmp = line.split()
onset_x = float(tmp[0])
onset_y = float(tmp[1])
onset_z = float(tmp[2])
# xyz file
with open(working_dir + "/mol/" + xyz, "r") as f:
lines = f.readlines()
num_atoms = int(lines[0].strip())
atoms = [line.strip() for line in lines[2:2+num_atoms]]
sys_dict = {
"molecule": "inline",
"geometry": "\n".join(atoms), # XYZ format as a string
"basis_file": working_dir + "/cap_data/basis/" + input_basis,
"bohr_coordinates": unit == 'Bohr'
}
cap_system = pyopencap.System(sys_dict)
print("opencap")
print(cap_system.get_overlap_mat("pyscf"))
print("pyscf")
print(ovlp)
if not(cap_system.check_overlap_mat(ovlp, "pyscf")):
raise Exception(
"Provided cap basis does not match to the pyscf basis.")
cap_dict = {"cap_type": "box",
"cap_x": onset_x,
"cap_y": onset_y,
"cap_z": onset_z,
"Radial_precision": "16",
"angular_points": "590",
"thresh": 10}
pc = pyopencap.CAP(cap_system, cap_dict, norb)
cap_ao = pc.get_ao_cap(ordering="pyscf")
def write_matrix_to_file(matrix, size, file, cutoff=1e-15):
f = open(file, 'w')
@ -143,7 +187,7 @@ def write_matrix_to_file(matrix, size, file, cutoff=1e-15):
# Write all 1 electron quantities in files
# Ov,Nuc,Kin,x,y,z
# Ov,Nuc,Kin,x,y,z,CAP
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'])
@ -156,8 +200,8 @@ 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'])
subprocess.call(['rm', '-f', working_dir + '/int/CAP.dat'])
write_matrix_to_file(cap_ao, norb, working_dir+'/int/CAP.dat')
def write_tensor_to_file(tensor, size, file_name, cutoff=1e-15):

107
basis_cap/install.py
View File

@ -1,107 +0,0 @@
#!/usr/bin/env python3
"""Installs the default basis set files from the BSE command-line tool in the required format for """
import basis_set_exchange as bse
conversion = {
"3-21g": "3-21G",
"4-31g": "4-31G",
"6-311g": "6-311G",
"6-311++g_2d_2p": "6-311++G(2d,2p)",
"6-311g_2df_2pd": "6-311G(2df,2pd)",
"6-311++g_3df_3pd": "6-311++G(3df,3pd)",
"6-311+g_star": "6-311+G*",
"6-311G_star": "6-311G*",
"6-311++g_star_star": "6-311++G**",
"6-311G_star_star": "6-311G**",
"6-31g": "6-31G",
"6-31+g": "6-31+G",
"6-31++g": "6-31++G",
"6-31g_3df_3pd": "6-31G(3df,3pd)",
"6-31++g_star": "6-31++G*",
"6-31+g_star": "6-31+G*",
"6-31g_star": "6-31G*",
"6-31++g_star_star": "6-31++G**",
"6-31g_star_star": "6-31G**",
"ano2_ames": "NASA Ames ANO2",
"ano_ames": "NASA Ames ANO",
"ano-rcc": "ANO-RCC",
"apr-cc-pv_q+d_z": "apr-cc-pV(Q+d)Z",
"aug-cc-pcv5z": "aug-cc-pCV5Z",
"aug-cc-pcvdz": "aug-cc-pCVDZ",
"aug-cc-pcvqz": "aug-cc-pCVQZ",
"aug-cc-pcvtz": "aug-cc-pCVTZ",
"aug-cc-pv5z": "aug-cc-pV5Z",
"aug-cc-pv6z": "aug-cc-pV6Z",
"aug-cc-pvdz": "aug-cc-pVDZ",
"aug-cc-pvqz": "aug-cc-pVQZ",
"aug-cc-pvtz": "aug-cc-pVTZ",
"aug-cc-pwcv5z": "aug-cc-pwCV5Z",
"aug-cc-pwcvdz": "aug-cc-pwCVDZ",
"aug-cc-pwcvqz": "aug-cc-pwCVQZ",
"aug-cc-pwcvtz": "aug-cc-pwCVTZ",
"cc-pcv5z": "cc-pCV5Z",
"cc-pcvdz": "cc-pCVDZ",
"cc-pcvqz": "cc-pCVQZ",
"cc-pcvtz": "cc-pCVTZ",
"cc-pv5z": "cc-pV5Z",
"cc-pv6z": "cc-pV6Z",
"cc-pv8z": "cc-pV8Z",
"cc-pv9z": "cc-pV9Z",
"cc-pvdz": "cc-pVDZ",
"cc-pvqz": "cc-pVQZ",
"cc-pvtz": "cc-pVTZ",
"cc-pwcv5z": "cc-pwCV5Z",
"cc-pwcvdz": "cc-pwCVDZ",
"cc-pwcvqz": "cc-pwCVQZ",
"cc-pwcvtz": "cc-pwCVTZ",
"d-aug-cc-pv5z": "d-aug-cc-pV5Z",
"d-aug-cc-pv6z": "d-aug-cc-pV6Z",
"d-aug-cc-pvdz": "d-aug-cc-pVDZ",
"d-aug-cc-pvqz": "d-aug-cc-pVQZ",
"d-aug-cc-pvtz": "d-aug-cc-pVTZ",
"def2-qzvp": "def2-QZVP",
"def2-qzvpd": "def2-QZVPD",
"def2-qzvpp": "def2-QZVPP",
"def2-qzvppd": "def2-QZVPPD",
"def2-svp": "def2-SVP",
"def2-sv_p": "def2-SV(P)",
"def2-svpd": "def2-SVPD",
"def2-tzvp": "def2-TZVP",
"def2-tzvpd": "def2-TZVPD",
"def2-tzvpp": "def2-TZVPP",
"def2-tzvppd": "def2-TZVPPD",
"jul-cc-pv_d+d_z": "jul-cc-pV(D+d)Z",
"jul-cc-pv_q+d_z": "jul-cc-pV(Q+d)Z",
"jul-cc-pv_t+d_z": "jul-cc-pV(T+d)Z",
"jun-cc-pv_d+d_z": "jun-cc-pV(D+d)Z",
"jun-cc-pv_q+d_z": "jun-cc-pV(Q+d)Z",
"jun-cc-pv_t+d_z": "jun-cc-pV(T+d)Z",
"lanl08": "LANL08",
"lanl08+": "LANL08+",
"lanl08d": "LANL08(d)",
"lanl08f": "LANL08(f)",
"m6-31g": "m6-31G",
"may-cc-pv_q+d_z": "may-cc-pV(Q+d)Z",
"may-cc-pv_t+d_z": "may-cc-pV(T+d)Z",
"midi_bang": "MIDI!",
"midi_huzinaga": "MIDI",
"mini": "MINI",
"pv6z": "pV6Z",
"pv7z": "pV7Z",
"pvdz_ahlrichs": "Ahlrichs pVDZ",
"sto-2g": "STO-2G",
"sto-3g": "STO-3G",
"sto-4g": "STO-4G",
"sto-5g": "STO-5G",
"sto-6g": "STO-6G",
"sto-3g_star": "STO-3G*",
}
for filename in conversion.keys():
print(filename)
data = bse.get_basis(
conversion[filename], fmt='psi4', uncontract_general=True, uncontract_spdf=True)
with open(filename, 'w') as f:
f.write(data)

1
cap_onsets/CO
View File

@ -1 +0,0 @@
2.76 2.76 4.97

1
cap_onsets/N2
View File

@ -1 +0,0 @@
2.76 2.76 4.88

2
crhf_test/methods.test
View File

@ -13,7 +13,7 @@
# G0F2 evGF2 qsGF2 ufGF2 G0F3 evGF3
F F F F F F
# G0W0 evGW qsGW ufG0W0 ufGW
F F F F F
T F F F F
# G0T0pp evGTpp qsGTpp ufG0T0pp
F F F F
# G0T0eh evGTeh qsGTeh

4
crhf_test/options.test
View File

@ -1,5 +1,7 @@
# HF: maxSCF thresh DIIS guess mix shift stab search
1000 0.00001 5 1 0.0 0.0 F F
# CAP: eta
0.0015
# MP: reg
F
# CC: maxSCF thresh DIIS
@ -9,7 +11,7 @@
# GF: maxSCF thresh DIIS lin eta renorm reg
256 0.00001 5 F 0.0 0 F
# GW: maxSCF thresh DIIS lin eta TDA_W reg
256 0.00001 5 T 0.0 F F
256 0.00001 5 F 0.0 F F
# GT: maxSCF thresh DIIS lin eta TDA_T reg
256 0.00001 5 F 0.0 F F
# ACFDT: AC Kx XBS

2
crhf_test/run_test.sh
View File

@ -5,4 +5,4 @@ cp ./options.test ../input/options
basis=$2
molecule=$1
cd ..
python3 PyDuck.py -x $molecule -b $basis -c 0 -m 1
python3 PyDuck.py -x $molecule -b $basis -c 0 -m 1 --use_cap

100
int/cap/CO/sto-3g/CAP.dat
View File

@ -1,100 +0,0 @@
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

100
int/cap/N2/sto-3g/CAP.dat
View File

@ -1,100 +0,0 @@
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

183
src/GW/cRG0W0.f90
View File

@ -1,5 +1,5 @@
subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,CAP,dipole_int,eHF)
! Perform G0W0 calculation
@ -37,6 +37,7 @@ subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TD
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: CAP(nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
double precision,intent(in) :: eHF(nOrb)
@ -46,24 +47,28 @@ subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TD
logical :: plot_self = .false.
logical :: dRPA_W
integer :: isp_W
integer :: p
double precision :: flow
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcGM
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: Re_SigC(:)
double precision,allocatable :: Im_SigC(:)
double precision,allocatable :: Re_Z(:)
double precision,allocatable :: Im_Z(:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:)
double precision,allocatable :: eGWlin(:)
double precision,allocatable :: eGW(:)
! Output variables
double precision,allocatable :: Re_eGWlin(:)
double precision, allocatable :: Im_eGWlin(:)
double precision,allocatable :: Re_eGW(:)
double precision,allocatable :: Im_eGW(:)
double precision, allocatable :: e_cap(:)
! Hello world
@ -88,17 +93,20 @@ subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TD
flow = 500d0
if(doSRG) then
! Not implemented
write(*,*) '*** SRG regularized G0W0 scheme ***'
write(*,*) '!!! No SRG with cRG0W0 !!!'
write(*,*)
end if
! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),SigC(nOrb),Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), &
eGW(nOrb),eGWlin(nOrb))
allocate(Aph(nS,nS),Bph(nS,nS),Re_SigC(nOrb),Im_SigC(nOrb),Re_Z(nOrb),Im_Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), &
Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),e_cap(nOrb))
do p = 1, nOrb
e_cap(p) = CAP(p,p)
end do
!-------------------!
! Compute screening !
!-------------------!
@ -119,37 +127,33 @@ subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TD
!------------------------!
! Compute GW self-energy !
!------------------------!
!!! HERE I HAVE TO WORK !!!
if(doSRG) then
write(*,*) "SRG for this method is not yet implemented !"
!call RGW_SRG_self_energy_diag(flow,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
else
call cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
end if
call cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z,CAP)
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
! Linearized or graphical solution?
eGWlin(:) = eHF(:) + Z(:)*SigC(:)
Re_eGWlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:)
Im_eGWlin(:) = Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:)
do p = 1, nOrb
Im_eGWlin(p) = Im_eGWlin(p) + CAP(p,p)
end do
if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
eGW(:) = eGWlin(:)
Re_eGW(:) = Re_eGWlin(:)
Im_eGW(:) = Im_eGWlin(:)
else
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call RGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z)
call cRGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,eHF,e_cap,Om,rho,Re_eGWlin,Im_eGWlin,eHF,e_cap,Re_eGW,Im_eGW,Re_Z,Im_Z)
end if
! Plot self-energy, renormalization factor, and spectral function
@ -162,7 +166,7 @@ subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TD
! Compute the RPA correlation energy
call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,Re_eGW,ERI,Aph)
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -171,72 +175,71 @@ subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TD
! Dump results !
!--------------!
call print_RG0W0(nOrb,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
call print_cRG0W0(nOrb,nO,eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM,CAP)
!---------------------------!
! Perform phBSE calculation !
!---------------------------!
if(dophBSE) then
call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, &
nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the BSE correlation energy via the adiabatic connection fluctuation dissipation theorem
if(doACFDT) then
call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
end if
!---------------------------!
! Perform ppBSE calculation !
!---------------------------!
if(doppBSE) then
call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
! Testing zone
if(dotest) then
call dump_test_value('R','G0W0 correlation energy',EcRPA)
call dump_test_value('R','G0W0 HOMO energy',eGW(nO))
call dump_test_value('R','G0W0 LUMO energy',eGW(nO+1))
end if
!
! if(dophBSE) then
!
! call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, &
! nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE)
!
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! ! Compute the BSE correlation energy via the adiabatic connection fluctuation dissipation theorem
!
! if(doACFDT) then
!
! call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eHF,Re_eGW,EcBSE)
!
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! end if
!
! end if
!
!!---------------------------!
!! Perform ppBSE calculation !
!!---------------------------!
!
! if(doppBSE) then
!
! call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE)
!
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! end if
!
!! Testing zone
!
! if(dotest) then
!
! call dump_test_value('R','G0W0 correlation energy',EcRPA)
! call dump_test_value('R','G0W0 HOMO energy',Re_eGW(nO))
! call dump_test_value('R','G0W0 LUMO energy',Re_eGW(nO+1))
!
! end if
!
end subroutine

91
src/GW/cRGW_self_energy.f90
View File

@ -1,91 +0,0 @@
subroutine RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
! Compute diagonal of the correlation part of the self-energy and the renormalization factor
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: i,a,p,m
double precision :: num,eps
! Output variables
double precision,intent(out) :: Sig(nBas)
double precision,intent(out) :: Z(nBas)
double precision,intent(out) :: EcGM
! Initialize
Sig(:) = 0d0
Z(:) = 0d0
!----------------!
! GW self-energy !
!----------------!
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do m=1,nS
eps = e(p) - e(i) + Om(m)
num = 2d0*rho(p,i,m)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do m=1,nS
eps = e(p) - e(a) - Om(m)
num = 2d0*rho(p,a,m)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
! Galitskii-Migdal correlation energy
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do m=1,nS
eps = e(a) - e(i) + Om(m)
num = 4d0*rho(a,i,m)**2
EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do
end do
end do
! Compute renormalization factor from derivative
Z(:) = 1d0/(1d0 - Z(:))
end subroutine

12
src/HF/cRHF.f90
View File

@ -1,5 +1,5 @@
subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
nBas,nO,S,T,V,ERI,X,ERHF,eHF,c,P,F)
nBas,nO,S,T,V,ERI,CAP,X,ERHF,eHF,c,P,F)
! Perform complex restricted Hartree-Fock calculation
@ -23,6 +23,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
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) :: CAP(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
@ -38,8 +39,6 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
double precision :: rcond
double precision,external :: trace_matrix
double precision :: eta
double precision,allocatable :: W(:,:)
complex*16,allocatable :: J(:,:)
complex*16,allocatable :: K(:,:)
complex*16,allocatable :: cp(:,:)
@ -69,25 +68,20 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
! Useful quantities
nBasSq = nBas*nBas
eta = 0.01d0
! Memory allocation
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(:,:)
! Define core Hamiltonian with CAP part
Hc(:,:) = cmplx(T+V,W,kind=8)
Hc(:,:) = cmplx(T+V,CAP,kind=8)
! Guess coefficients and density matrix
call complex_mo_guess(nBas,nBas,guess_type,S,Hc,X,c)

15
src/QuAcK/QuAcK.f90
View File

@ -13,6 +13,7 @@ program QuAcK
logical :: dophRPA,dophRPAx,docrRPA,doppRPA
logical :: doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3
logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
logical :: docG0W0
logical :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh
integer :: nNuc
@ -31,6 +32,7 @@ program QuAcK
double precision,allocatable :: V(:,:)
double precision,allocatable :: Hc(:,:)
double precision,allocatable :: X(:,:),X_tmp(:,:)
double precision, allocatable :: CAP(:,:)
double precision,allocatable :: dipole_int_AO(:,:,:)
double precision,allocatable :: Uvec(:,:), Uval(:)
@ -41,6 +43,8 @@ program QuAcK
double precision :: thresh_HF,level_shift,mix
integer :: guess_type
double precision :: eta_cap
logical :: reg_MP
logical :: switch_hpc
@ -123,6 +127,7 @@ program QuAcK
doG0W0,doevGW,doqsGW,doufG0W0,doufGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, &
doG0T0eh,doevGTeh,doqsGTeh, &
docG0W0, &
doRtest,doUtest,doGtest)
!--------------------------!
@ -131,6 +136,7 @@ program QuAcK
call read_options(working_dir, &
maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, &
eta_cap, &
reg_MP, &
maxSCF_CC,thresh_CC,max_diis_CC, &
TDA,spin_conserved,spin_flip, &
@ -182,12 +188,14 @@ program QuAcK
allocate(V(nBas,nBas))
allocate(Hc(nBas,nBas))
allocate(dipole_int_AO(nBas,nBas,ncart))
allocate(CAP(nBas,nBas))
! Read integrals
call wall_time(start_int)
call read_1e_integrals(working_dir,nBas,S,T,V,Hc)
if (docRHF .or. docG0W0) call read_CAP_integrals(nBas,CAP) ! Add different cases if needed
CAP(:,:) = -eta_cap * CAP(:,:)
call read_dipole_integrals(working_dir,nBas,dipole_int_AO)
call wall_time(end_int)
@ -250,8 +258,9 @@ program QuAcK
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
docG0W0, &
nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
S,T,V,Hc,CAP,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
@ -327,5 +336,5 @@ program QuAcK
if (allocated(V)) deallocate(V)
if (allocated(Hc)) deallocate(Hc)
if (allocated(dipole_int_AO)) deallocate(dipole_int_AO)
!if (allocated(S)) deallocate(S)
if (allocated(S)) deallocate(S)
end program

27
src/QuAcK/RQuAcK.f90
View File

@ -3,8 +3,9 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
docG0W0, &
nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
S,T,V,Hc,CAP_AO,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
@ -33,6 +34,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
logical,intent(in) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp
logical,intent(in) :: doG0T0eh,doevGTeh,doqsGTeh
logical,intent(in) :: docG0W0
integer,intent(in) :: nNuc,nBas,nOrb
integer,intent(in) :: nC
@ -47,6 +49,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: CAP_AO(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
@ -107,6 +110,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
complex*16,allocatable :: complex_FHF(:,:)
double precision :: ERHF
complex*16 :: complex_ERHF
double precision,allocatable :: CAP_MO(:,:)
double precision,allocatable :: dipole_int_MO(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
@ -131,6 +135,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
allocate(complex_FHF(nBas,nBas))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
allocate(FHF(nBas,nBas))
allocate(CAP_MO(nOrb,nOrb))
allocate(dipole_int_MO(nOrb,nOrb,ncart))
allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
call wall_time(start_int)
@ -174,7 +179,7 @@ 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,ENuc, &
nBas,nO,S,T,V,ERI_AO,X,complex_ERHF,complex_eHF,complex_cHF,complex_PHF,complex_FHF)
nBas,nO,S,T,V,ERI_AO,CAP_AO,X,complex_ERHF,complex_eHF,complex_cHF,complex_PHF,complex_FHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
@ -202,7 +207,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
! 4-index transform
call AOtoMO_ERI_RHF(nBas,nOrb,cHF,ERI_AO,ERI_MO)
if (docG0W0) call AOtoMO(nBas,nOrb,cHF,CAP_AO,CAP_MO) ! Add the different cases for cGW methods when implemented
call wall_time(end_AOtoMO)
t_AOtoMO = end_AOtoMO - start_AOtoMO
@ -343,15 +348,14 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
! GW module !
!-----------!
doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW
doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW .or. docG0W0
if(doGW) then
call wall_time(start_GW)
call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF_GW,thresh_GW,max_diis_GW, &
call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF_GW,thresh_GW,max_diis_GW, &
doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T, &
V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -360,6 +364,13 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
end if
!------------!
! cGW module !
!------------!
! IMPLEMENT LATER TO TREAT EVERYTHING COMPLEX, i.e. start from complex HF orbitals
!-----------------!
! T-matrix module !
!-----------------!
@ -391,7 +402,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
deallocate(dipole_int_MO)
deallocate(ERI_MO)
deallocate(ERI_AO)
deallocate(CAP_MO)
deallocate(complex_eHF)
deallocate(complex_cHF)
deallocate(complex_PHF)

11
src/QuAcK/read_methods.f90
View File

@ -10,6 +10,7 @@ subroutine read_methods(working_dir, &
doG0W0,doevGW,doqsGW,doufG0W0,doufGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, &
doG0T0eh,doevGTeh,doqsGTeh, &
docG0W0, &
doRtest,doUtest,doGtest)
! Read desired methods
@ -32,6 +33,7 @@ subroutine read_methods(working_dir, &
logical,intent(out) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
logical,intent(out) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp
logical,intent(out) :: doG0T0eh,doevGTeh,doqsGTeh
logical,intent(out) :: docG0W0
logical,intent(out) :: doRtest,doUtest,doGtest
@ -175,6 +177,8 @@ subroutine read_methods(working_dir, &
if(ans4 == 'T') doufG0W0 = .true.
if(ans5 == 'T') doufGW = .true.
! Read GTpp methods
doG0T0pp = .false.
@ -201,6 +205,13 @@ subroutine read_methods(working_dir, &
if(ans2 == 'T') doevGTeh = .true.
if(ans3 == 'T') doqsGTeh = .true.
! Read CAP-GW methods
docG0W0 = .false.
read(1,*)
read(1,*) ans1
if(ans1 == 'T') docG0W0 = .true.
! Read test
doRtest = .false.

8
src/QuAcK/read_options.f90
View File

@ -1,5 +1,6 @@
subroutine read_options(working_dir, &
maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, &
eta_cap, &
reg_MP, &
maxSCF_CC,thresh_CC,max_diis_CC, &
TDA,spin_conserved,spin_flip, &
@ -29,6 +30,8 @@ subroutine read_options(working_dir,
logical,intent(out) :: dostab
logical,intent(out) :: dosearch
double precision,intent(out) :: eta_cap
logical,intent(out) :: reg_MP
integer,intent(out) :: maxSCF_CC
@ -113,6 +116,11 @@ subroutine read_options(working_dir,
if(ans1 == 'T') dostab = .true.
if(ans2 == 'T') dosearch = .true.
! Read cap options
read(1,*)
read(1,*) eta_cap
! Read MPn options
reg_MP = .false.

8
src/utils/read_CAP_integrals.f90
View File

@ -13,7 +13,7 @@ subroutine read_CAP_integrals(nBas,W)
logical :: debug
integer :: mu,nu
double precision :: wx,wy,wz
double precision :: wxyz
! Output variables
@ -28,9 +28,9 @@ subroutine read_CAP_integrals(nBas,W)
W(:,:) = 0d0
do
read(31,*,end=31) mu,nu,wx,wy,wz
W(mu,nu) = wx + wy + wz
W(nu,mu) = wx + wy + wz
read(31,*,end=31) mu,nu,wxyz
W(mu,nu) = wxyz
W(nu,mu) = wxyz
end do
31 close(unit=31)