Skip to content

GitLab

Commit 147d4c16 authored by Niccolo.Ricardi's avatar Niccolo.Ricardi
Browse files

FnT_general starts from last iteration

parent 257d2665
Hide whitespace changes
Inline Side-by-side
FDET/freeze_and_thaw.py
View file @ 147d4c16
......@@ -26,13 +26,13 @@ class NotConvergedError(Exception):
pass
def print_iteration(it, lw=80, opt="FT"):
dict_={"FT":"Freeze-And-Thaw", "MC": "Conventional FDET"}
dict_={"FT":"Freeze-and-Thaw", "MC": "Conventional FDET"}
if opt not in dict_.keys():
dict_ [opt]="Unknown kind of cycle"
sw = 20
q,m = divmod(lw-sw,2)
print("*"*lw)
print(" "*q, "{0} ##{1}".format(dict_[opt],it))
print(" "*q, "{0} #{1}".format(dict_[opt],it))
print("*"*lw)
def print_banner(string, lw=80):
......@@ -293,26 +293,27 @@ def get_embedded_energy(outfile, with_ccp=False):
# nbas["nbasB"] = parsed[index[2]][1]
# return nbas
def get_nbas(file, with_ccp=False, index=[2,3,4]): #json not always named as outfile. done to work with "file,out" and "file.json"
def get_nbas(filepath, with_ccp=False, index=[2,3,4]): #json not always named as outfile. done to work with "file,out" and "file.json"
""" Get number of basis functions"""
nbas = {"nbasA" : 0, "nbasB": 0, "nbasAB": 0}
dir_,file = os.path.split(filepath)[0], os.path.split(filepath)[1]
if with_ccp:
if file[-4:]=="json":
json_file=file
elif os.path.exists("ccp_{0}.json".format(file[:-4])):
json_file="ccp_{0}.json".format(file[:-4])
if filepath[-4:] == "json":
json_file = filepath
elif os.path.exists(os.path.join(dir_,"ccp_{0}.json".format(file[:-4]))):
json_file=os.path.join(dir_,"ccp_{0}.json".format(file[:-4]))
else:
import glob as g
json_files=g.glob("*.json")
if len(json_files)==1:
json_files = g.glob(os.path.join(dir_,"*.json"))
if len(json_files) == 1:
json_file=json_files[0]
else:
out = ccp.Parser(file, to_console=False,
out = ccp.Parser(filepath, to_console=False,
to_file=True, to_json=True,
log_file="ccp_{0}.log".format(file[:-4]),
json_file="ccp_{0}.json".format(file[:-4]))
log_file=os.path.join(dir_,"ccp_{0}.log".format(file[:-4])),
json_file=os.path.join(dir_,"ccp_{0}.json".format(file[:-4])))
json_file="ccp_{0}.json".format(file[:-4])
json_file=os.path.join(dir_,"ccp_{0}.json".format(file[:-4]))
parsed = load_json(json_file)
nbas["nbasAB"] = parsed["nbas"][index[0]][0]
nbas["nbasA"] = parsed["nbas"][index[1]][0]
......@@ -742,89 +743,74 @@ def freezeAndThaw_general(specs, queue, specs_B=None, queue_B=None, method_A="HF
#------------------------------------------inp--------------------------------
# FIRST JOB
#--------------------------------------------------------------------------
print_iteration(it)
# Prepare first input
t_idx = 0 if not B_init else 1
curr_inp = iterFilName+".in"
curr_job = "FT-{0}".format(it)
# copy density matrix (B), if provided
if B_init:
shutil.copy(os.path.join(root, densMat_B), path_cur)
try:
E_embA = get_embedded_energy(iterFilName+".out", with_ccp=has_ccp)
energy.append(E_embA)
print_iteration(it)
print("It was already done!!Yay!")
except:
# Prepare first input
t_idx = 0 if not B_init else 1
curr_inp = iterFilName+".in"
curr_job = "FT-{0}".format(it)
templates = get_templates(method_A, isXC)
# copy density matrix (B), if provided
if B_init:
shutil.copy(os.path.join(root, densMat_B), path_cur)
templates = get_templates(method_A, isXC)
if method_B == "MP2":
if isXC:
from CCJob.Templates import ADCin_MP2
templates[0]=ADCin_MP2
else:
from CCJob.Templates import ADCin_MP2_X_C
templates[0]=ADCin_MP2_X_C
if specs["expansion"]=="ME":
specs["frag_b_bse"]=specs["frag_b"]
elif specs["expansion"]=="SE":
specs["frag_b_bse"] = "@"+"\n@".join(specs["frag_a"].split("\n"))+"\n"+specs["frag_b"]
if method_B == "MP2":
if isXC:
from CCJob.Templates import ADCin_MP2
templates[0]=ADCin_MP2
else:
from CCJob.Templates import ADCin_MP2_X_C
templates[0]=ADCin_MP2_X_C
if specs["expansion"]=="ME":
specs["frag_b_bse"]=specs["frag_b"]
elif specs["expansion"]=="SE":
specs["frag_b_bse"] = "@"+"\n@".join(specs["frag_a"].split("\n"))+"\n"+specs["frag_b"]
inp = ccj.Input.from_template(templates[t_idx], curr_inp, **specs)
# Prepare job
queue["jobname"] = curr_job
job = ccj.Job(inp, **queue)
if q_custom != None:
job.set_custom_options(**q_custom)
# Run the first job
job.run()
# get embedded energy of active subsystem
E_embA = get_embedded_energy(iterFilName+".out", with_ccp=has_ccp)
energy.append(E_embA)
inp = ccj.Input.from_template(templates[t_idx], curr_inp, **specs)
# --- RE-ORDER DENSITY MATRICES ---
A = read_density(densMat_from_FDE, debug=True)
hd, alp = method_B=="HF" and not B_init, method_B=="HF" and not B_init
B = read_density(densMat_B, header=hd, alpha_only=alp)
# transform each nAO x nAO matrix if supermolecular expansion
if specs["expansion"] == "SE":
# get number of basis functions from the basis initializer block (constructor)
indx=[3,4,5] if specs["expansion"]=="SE" and method_A=="MP2" and not B_init else [2,3,4]
nbas = get_nbas(iterFilName+".out", with_ccp=has_ccp,index=indx)
nbasAB = nbas["nbasAB"]
nbasA = nbas["nbasA"]
nbasB = nbas["nbasB"]
print("""-- Found number of basis functions:
.nbasAB = {0}
.nbasA = {1}
.nbasB = {2}
""".format(nbasAB, nbasA, nbasB))
A_tf = np.zeros((2, nbasAB, nbasAB)) # density that serves as rhoB in next it
B_tf = np.zeros((2, nbasAB, nbasAB))
# Prepare job
queue["jobname"] = curr_job
job = ccj.Job(inp, **queue)
if q_custom != None:
job.set_custom_options(**q_custom)
# Run the first job
job.run()
# get embedded energy of active subsystem
E_embA = get_embedded_energy(iterFilName+".out", with_ccp=has_ccp)
energy.append(E_embA)
# make a copy of A's specs dictionary
specsA = dict(**specs)
# get number of basis functions from the basis initializer block (constructor)
indx=[3,4,5] if specs["expansion"]=="SE" and method_A=="MP2" and not B_init else [2,3,4]
nbas = get_nbas(iterFilName+".out", with_ccp=has_ccp,index=indx)
nbasAB = nbas["nbasAB"]
nbasA = nbas["nbasA"]
nbasB = nbas["nbasB"]
print("""-- Found number of basis functions:
.nbasAB = {0}
.nbasA = {1}
.nbasB = {2}
""".format(nbasAB, nbasA, nbasB))
if specs["expansion"]=="SE":
number_line = [i for i in range(nbasAB)]
order_to_B = number_line[nbasA:] + number_line[:nbasA]
order_to_A = number_line[nbasB:] + number_line[:nbasB]
for i in range(2):
A_tf[i] = transform_AOmatrix(A[i], order=order_to_B)
B_tf[i] = transform_AOmatrix(B[i], order=order_to_B)
else:
A_tf = A
B_tf = B
order_to_A = number_line[nbasB:] + number_line[:nbasB]
# adjust iteration counter
it += 1
# make a copy of A's specs dictionary
specsA = dict(**specs)
it+=1
#------------------------------------------------------------------------------
# ALL OTHER JOBS
#------------------------------------------------------------------------------
difference_to_last = 1
while abs(difference_to_last) > thresh_fnt and it < max_iter:
# iteration variables
is_A = bool((it+1) % 2)
was_A = bool((it) % 2)
print_iteration(it)
curr_inp = iterFilName+".in"
curr_job = "FT-{0}".format(it)
......@@ -834,73 +820,82 @@ def freezeAndThaw_general(specs, queue, specs_B=None, queue_B=None, method_A="HF
if not os.path.exists(path_cur):
os.makedirs(path_cur)
os.chdir(path_cur)
# save transformed density matrix from last iteration
if is_A:
save_density("Densmat_B.txt", B_tf)
save_density("Densmat_A.txt", A_tf)
else:
save_density("Densmat_B.txt", A_tf)
save_density("Densmat_A.txt", B_tf)
# reverse order (embedded species is always "frag_a")
if is_A:
specs.update(specsA)
else:
specs["frag_a"] = specsA["frag_b"]
specs["frag_b"] = specsA["frag_a"]
if hasElConf:
specs["charge_a"] = specsA["charge_b"]
specs["charge_b"] = specsA["charge_a"]
specs["multiplicity_a"] = specsA["multiplicity_b"]
specs["multiplicity_b"] = specsA["multiplicity_a"]
# Prepare input & job
if not is_A and specs_B != None and queue_B != None:
queue_B["jobname"] = curr_job
inp = ccj.Input.from_template(templates[2], curr_inp,
frag_a=specsA["frag_b"],
frag_b=specsA["frag_a"], **specs_B)
job = ccj.Job(inp, **queue_B)
else:
queue["jobname"] = curr_job
inp = ccj.Input.from_template(templates[2], curr_inp, **specs)
job = ccj.Job(inp, **queue)
if q_custom != None:
job.set_custom_options(**q_custom)
# Run the first job
job.run()
# get embedded energy of active subsystem
E_embA = get_embedded_energy(iterFilName+".out", with_ccp=has_ccp)
energy.append(E_embA)
try:
E_embA = get_embedded_energy(iterFilName+".out", with_ccp=has_ccp)
energy.append(E_embA)
print_iteration(it)
print("It was already done!!Yay!")
it+=1
except:
# get number of basis functions from the basis initializer block (constructor)
#for ME the number of basis is not used. Kept as check
templates = get_templates(method_A, isXC)
prev_it=os.path.join("..",iterDirName+str(it-1))
# --- RE-ORDER DENSITY MATRICES ---
densMat_B = "Densmat_B.txt" if it>1 else densMat_B
A = read_density(os.path.join(prev_it,densMat_from_FDE), debug=True) #what is called A in the previous iteration
hd, alp = (not B_init), (not B_init)
hd,alp = False if it>1 else hd, False if it>1 else alp
B = read_density(os.path.join(prev_it,densMat_B), header=hd, alpha_only=alp) #what is called B in the previous iteration
# transform each nAO x nAO matrix if supermolecular expansion
if specs["expansion"] == "SE":
A_tf = np.zeros((2, nbasAB, nbasAB))
B_tf = np.zeros((2, nbasAB, nbasAB))
for i in range(2):
order= order_to_B if was_A else order_to_A
A_tf[i] = transform_AOmatrix(A[i], order=order)
B_tf[i] = transform_AOmatrix(B[i], order=order)
else:
A_tf = A
B_tf = B
# print current residual
if is_A:
difference_to_last = energy[it] - energy[it-2]
print("-- Current residual: {0: .8e}".format(difference_to_last))
# read embedded density from file: D.shape = (2, nAO, nAO)
A = read_density(densMat_from_FDE)
if specs["expansion"] == "SE":
for i in range(2):
if is_A:
A_tf[i] = transform_AOmatrix(A[i], order=order_to_B)
B_tf[i] = transform_AOmatrix(B_tf[i], order=order_to_B)
else:
A_tf[i] = transform_AOmatrix(A_tf[i], order=order_to_A)
B_tf[i] = transform_AOmatrix(A[i], order=order_to_A)
else:
if is_A:
A_tf = A
# save transformed density matrix from last iteration
save_density("Densmat_A.txt", B_tf)
save_density("Densmat_B.txt", A_tf)
# reverse order (embedded species is always "frag_a")
if not was_A:
try:
specs.update(specsA)
except:
specsA = dict(**specs)
else:
B_tf = A
# increase iteration counter
it += 1
# END WHILE LOOP
specs["frag_a"] = specsA["frag_b"]
specs["frag_b"] = specsA["frag_a"]
if hasElConf:
specs["charge_a"] = specsA["charge_b"]
specs["charge_b"] = specsA["charge_a"]
specs["multiplicity_a"] = specsA["multiplicity_b"]
specs["multiplicity_b"] = specsA["multiplicity_a"]
# Prepare input & job
if was_A and specs_B != None and queue_B != None:
queue_B["jobname"] = curr_job
inp = ccj.Input.from_template(templates[2], curr_inp,
frag_a=specsA["frag_b"],
frag_b=specsA["frag_a"], **specs_B)
job = ccj.Job(inp, **queue_B)
else:
queue["jobname"] = curr_job
inp = ccj.Input.from_template(templates[2], curr_inp, **specs)
job = ccj.Job(inp, **queue)
if q_custom != None:
job.set_custom_options(**q_custom)
# Run the first job
job.run()
# get embedded energy of active subsystem
E_embA = get_embedded_energy(iterFilName+".out", with_ccp=has_ccp)
energy.append(E_embA)
# print current residual
if not was_A:
difference_to_last = energy[it] - energy[it-2]
print("-- Current residual: {0: .8e}".format(difference_to_last))
# increase iteration counter
it += 1
# END WHILE LOOP
if it >= max_iter:
raise NotConvergedError("Freeze-and-Thaw cycles did not converge in "
"{0} iterations!".format(max_iter))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment