Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
103 changes: 83 additions & 20 deletions ash/interfaces/interface_dlfind.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ def DLFIND_optimizer(jobtype=None, theory=None, fragment=None, fragment2=None, c
hessian_choice="numfreq", inithessian=0,
numfreq_npoint=1, numfreq_displacement=0.005, numfreq_hessatoms=None,
numfreq_force_projection=None, print_atoms_list=None,
force_noPBC=False, PBC_format_option='CIF'):
force_noPBC=False, PBC_format_option='CIF', dump=None, restart=None):
"""
Wrapper function around DLFIND_optimizerClass
"""
Expand All @@ -46,7 +46,7 @@ def DLFIND_optimizer(jobtype=None, theory=None, fragment=None, fragment2=None, c
hessian_choice=hessian_choice, inithessian=inithessian,
numfreq_npoint=numfreq_npoint,numfreq_displacement=numfreq_displacement,
numfreq_hessatoms=numfreq_hessatoms,numfreq_force_projection=numfreq_force_projection,
print_atoms_list=print_atoms_list,
print_atoms_list=print_atoms_list, dump=dump, restart=restart,
force_noPBC=force_noPBC, PBC_format_option=PBC_format_option)

# If NumGrad then we wrap theory object into NumGrad class object
Expand All @@ -72,7 +72,7 @@ def __init__(self,jobtype=None, fragment=None, fragment2=None, theory=None, char
hessian_choice='numfreq', inithessian=None,
numfreq_npoint=1,numfreq_displacement=0.005,numfreq_force_projection=None,
numfreq_hessatoms=None, print_atoms_list=None,
force_noPBC=False, PBC_format_option='CIF'):
force_noPBC=False, PBC_format_option='CIF', dump=None, restart=None):

print_line_with_mainheader("DLFIND_optimizer initialization")
print()
Expand Down Expand Up @@ -180,6 +180,11 @@ def __init__(self,jobtype=None, fragment=None, fragment2=None, theory=None, char
self.tolerance_e=tolerance_e
# NEB
self.nimage=nimage

# DL-FIND restart/checkpoint options
self.dump=dump
self.restart=restart

#Dimer
self.delta=delta

Expand Down Expand Up @@ -231,7 +236,10 @@ def ash_e_g_func(coordinates, iimage, kiter, theory):
#PRINTING ACTIVE GEOMETRY IN EACH GEOMETRIC ITERATION
self.fragment.write_xyzfile(xyzfilename="Fragment-currentgeo.xyz")
if self.printlevel >= 1:
print(f"Current geometry (Å) in step {self.dlfind_opt_cycles} (print_atoms_list region)")
if self.icoord >= 100 and self.icoord < 150:
print(f"Current geometry (Å), NEB cycle {self.dlfind_neb_cycles + 1}, image {iimage} (print_atoms_list region)")
else:
print(f"Current geometry (Å) in step {self.dlfind_opt_cycles} (print_atoms_list region)")
print("---------------------------------------------------")
print_coords_for_atoms(R_phys, self.elems_phys, self.print_atoms_list)
print("")
Expand Down Expand Up @@ -279,7 +287,10 @@ def ash_e_g_func(coordinates, iimage, kiter, theory):
#PRINTING ACTIVE GEOMETRY IN EACH GEOMETRIC ITERATION
self.fragment.write_xyzfile(xyzfilename="Fragment-currentgeo.xyz")
if self.printlevel >= 1:
print(f"Current geometry (Å) in step {self.dlfind_opt_cycles} (print_atoms_list region)")
if self.icoord >= 100 and self.icoord < 150:
print(f"Current geometry (Å), NEB cycle {self.dlfind_neb_cycles + 1}, image {iimage} (print_atoms_list region)")
else:
print(f"Current geometry (Å) in step {self.dlfind_opt_cycles} (print_atoms_list region)")
print("---------------------------------------------------")
print_coords_for_atoms(coordinates_ang, self.fragment.elems, self.print_atoms_list)
print("")
Expand Down Expand Up @@ -371,7 +382,7 @@ def hess_func(coords):
# Create function to store results from DL-FIND
#@dlf_put_coords_wrapper
def store_results(a,nvar,switch, energy, coordinates, iam):
print("Called store_results with switch:", switch)
print("Called store_results with switch:", switch, flush=True)
if switch > 0:
coords = as_array(coordinates, (nvar,)).reshape(-1, 3)
coordinates_ang = coords*0.5291772109303
Expand All @@ -380,23 +391,30 @@ def store_results(a,nvar,switch, energy, coordinates, iam):
# Write out NEB path if switch -1
if self.icoord >= 100 and self.icoord < 150 and switch == -1:
self.dlfind_neb_cycles+=1
print("="*70)
print(f"DLFIND NEB-OPTIMIZATION CYCLE {self.dlfind_neb_cycles}")
print("="*70)
print("="*70, flush=True)
print(f"DLFIND NEB-OPTIMIZATION CYCLE {self.dlfind_neb_cycles}", flush=True)
print("="*70, flush=True)

for imageid in list(range(1,self.nimage)):
if imageid not in self.NEB_energies_dict or imageid not in self.NEB_geometries:
print(f"Error: Missing ASH NEB data for image {imageid}")
ashexit()
#print("Switch -1, writing out current NEB-path without CI")
with open("DLFIND_NEBpath_all.xyz", "a") as trajfile:
for imageid in list(range(1,self.nimage)):
trajfile.write(str(self.fragment.numatoms) + "\n")
trajfile.write(f"Image {imageid}. Energy: {self.NEB_energies_dict[imageid]} \n")
for el, cord in zip(self.fragment.elems, self.NEB_geometries[imageid]):
trajfile.write(el + " " + str(cord[0]) + " " + str(cord[1]) + " " + str(cord[2]) + "\n")
with open("DLFIND_NEBpath_current.xyz", "w") as trajfile:
with open("DLFIND_NEBpath_current.xyz.tmp", "w") as trajfile:
for imageid in list(range(1,self.nimage)):
trajfile.write(str(self.fragment.numatoms) + "\n")
trajfile.write(f"Image {imageid}. Energy: {self.NEB_energies_dict[imageid]} \n")
for el, cord in zip(self.fragment.elems, self.NEB_geometries[imageid]):
trajfile.write(el + " " + str(cord[0]) + " " + str(cord[1]) + " " + str(cord[2]) + "\n")

os.replace("DLFIND_NEBpath_current.xyz.tmp", "DLFIND_NEBpath_current.xyz")

# Writing out CI once it has been spawned
if nimage in self.NEB_geometries:
print("Writing out current trajectory for climbing image as: DLFIND_CIgeo_traj.xyz ")
Expand Down Expand Up @@ -673,24 +691,29 @@ def prepare_run(self):
iopt=self.iopt, maxcycle=self.maxcycle,tolerance=self.tolerance,
tolerance_e=self.tolerance_e, inithessian=self.inithessian,
nframe=nframe, nz=nz,
ncons=self.numcons, delta=self.delta,
ncons=self.numcons, delta=self.delta, dump=self.dump, restart=self.restart,
spec=self.spec, printl=self.printlevel, nimage=self.nimage)

# Delete old traj file before beginning
# Delete old traj files before beginning, unless this is a restart
remove_files=['DLFIND_opt_traj.xyz','DLFIND_dimertraj_1.xyz', 'DLFIND_dimertraj_2.xyz','DLFIND_dimertraj_3.xyz','DLFIND_NEBpath_current.xyz', 'DLFIND_NEBpath_all.xyz', 'DLFIND_CIgeo_traj.xyz']
print_if_level(f"Removing possible old files: {remove_files}", self.printlevel,2)
for rfile in remove_files:
try:
os.remove(rfile)
print_if_level(f"removed {rfile} ", self.printlevel,2)
except FileNotFoundError:
#print(f"file {rfile} not found")
pass

if self.restart in (None,0):
print_if_level(f"Removing possible old files: {remove_files}", self.printlevel,2)
for rfile in remove_files:
try:
os.remove(rfile)
print_if_level(f"removed {rfile} ", self.printlevel,2)
except FileNotFoundError:
pass
else:
print_if_level("Restart requested. Keeping existing DL-FIND files.", self.printlevel,2)

print_if_level(f"\nArguments passed to DL-FIND:", self.printlevel,2)
print_if_level(f"icoord: {self.icoord}", self.printlevel,2)
print_if_level(f"iopt: {self.iopt}", self.printlevel,2)
print_if_level(f"maxcycle: {self.maxcycle}", self.printlevel,2)
print_if_level(f"dump: {self.dump}", self.printlevel,2)
print_if_level(f"restart: {self.restart}", self.printlevel,2)
print_if_level(f"spec ({len(self.spec)}): {self.spec}", self.printlevel,2)
if self.icoord == 120:
print_if_level(f"NEB nimage: {self.nimage}", self.printlevel,2)
Expand All @@ -707,6 +730,46 @@ def run(self, theory=None, fragment=None, fragment2=None, constraints=None, char
self.NEB_energies_dict={}
self.NEB_geometries={}

# Restore ASH NEB bookkeeping on restart
if self.restart not in (None,0) and 100 <= self.icoord < 150:
# Restore previous ASH NEB cycle count
if os.path.isfile("DLFIND_NEBpath_all.xyz"):
with open("DLFIND_NEBpath_all.xyz") as f:
self.dlfind_neb_cycles = sum(1 for line in f if line.startswith("Image 1."))
print(f"Restart: restored ASH NEB cycle counter to {self.dlfind_neb_cycles}")

# Restore energies/geometries of the most recent complete NEB path
if os.path.isfile("DLFIND_NEBpath_current.xyz"):
print("Restart: restoring ASH NEB energies/geometries from DLFIND_NEBpath_current.xyz")
with open("DLFIND_NEBpath_current.xyz") as f:
lines = f.readlines()

i = 0
while i < len(lines):
natoms = int(lines[i].strip())
header = lines[i + 1].split()
imageid = int(header[1].rstrip("."))
energy = float(header[3])

coords = []
for j in range(natoms):
parts = lines[i + 2 + j].split()
coords.append([float(parts[1]), float(parts[2]), float(parts[3])])

self.NEB_energies_dict[imageid] = energy
self.NEB_geometries[imageid] = np.array(coords)

i += natoms + 2

print("Restored NEB images:", sorted(self.NEB_energies_dict.keys()))

expected_images = set(range(1,self.nimage))
restored_images = set(self.NEB_energies_dict.keys())
missing_images = expected_images - restored_images
if len(missing_images) > 0:
print(f"Error: Missing NEB restart data for images: {sorted(missing_images)}")
ashexit()

# Runcounter
self.runcounter=0

Expand Down
Loading