Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

IRC bug fix and NEB update #207

Merged
merged 2 commits into from
Jul 19, 2024
Merged

IRC bug fix and NEB update #207

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
194f581
irc update
hjnpark Jun 19, 2024
d92ef1b
verbose option is added for irc
hjnpark Jun 19, 2024
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
44 changes: 28 additions & 16 deletions geometric/optimize.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -450,7 +450,9 @@ def guess_g(self, g, H, disp):

def IRC_step(self):
self.farConstraints = self.IC.haveConstraints() and self.IC.maxConstraintViolation(self.X) > 1e-1
logger.info("IRC sub-step 1: Finding the pivot point (q*_{k+1})\n")

if self.params.verbose:
logger.info("IRC sub-step 1: Finding the pivot point (q*_{k+1})\n")

# Need to take a step towards the pivot point
self.IC.clearCache()
Expand All @@ -462,7 +464,7 @@ def IRC_step(self):
# Vector to the pivot point
if self.Iteration == 0:
# If it's the very first step, pick the eigenvector of the imaginary frequency and pick the direction
logger.info('\nFirst, following the imaginary mode vector\n')
logger.info('First, following the imaginary mode vector\n')
if self.TSWavenum[1] < 0:
raise IRCError("There are more than one imaginary vibrational mode. Please optimize the structure and try again.\n")
elif self.TSWavenum[0] > 0:
Expand Down Expand Up @@ -504,11 +506,13 @@ def IRC_step(self):
# Calculating sqrt(mass) weighted Cartesian coordinate
MWGMat_sqrt_inv, MWGMat_sqrt = self.IC.GInverse_SVD(X_pivot, sqrt=True, invMW=True)
mwdx_1 = np.dot(MWGMat_sqrt_inv, dy_to_pivot)
logger.info("Half step dy = %.5f\n" %np.linalg.norm(dy_to_pivot))
logger.info("Half step mw-dx = %.5f Bohr*sqrt(amu)\n" %np.linalg.norm(mwdx_1))

if self.params.verbose:
logger.info("Half step dy = %.5f\n" %np.linalg.norm(dy_to_pivot))
logger.info("Half step mw-dx = %.5f Bohr*sqrt(amu)\n" %np.linalg.norm(mwdx_1))

# We are at the pivot point
logger.info('\nIRC sub-step 2: Finding the next point (q_{k+1})\n')
logger.info('\nIRC sub-step 2: Finding the next point (q_{k+1})\n')
v1 = v.copy()
irc_sub_iteration = 0
irc_reset_iteration = 0
Expand Down Expand Up @@ -571,11 +575,12 @@ def IRC_step(self):
irc_sub_iteration += 1
p_prime += dq_new

logger.info('Angle between v1 and v2: %2.f \n' % deg)
logger.info('Half step dy = %.5f \n' % np.linalg.norm(p_prime))
logger.info('Half step mw-dx = %.5f Bohr*sqrt(amu)\n\n' % half_mwdx)
if self.params.verbose:
logger.info('Angle between v1 and v2: %2.f \n' % deg)
logger.info('Half step dy = %.5f \n' % np.linalg.norm(p_prime))
logger.info('Half step mw-dx = %.5f Bohr*sqrt(amu)\n\n' % half_mwdx)
logger.info('=> Total step dy = %.5f \n' % dy_norm)
logger.info('=> Total step mw-dx = %.5f Bohr*sqrt(amu)\n' % mwdx)
logger.info('=> Total step mw-dx = %.5f Bohr*sqrt(amu)\n\n' % mwdx)
self.Iteration += 1
return dy

Expand Down Expand Up @@ -744,7 +749,7 @@ def evaluate_IRC_step(self, params, step_state, criteria_met, IRC_converged):
return True, step_state

self.IC_check = False
if step_state in (StepState.Reject, StepState.Poor) and not self.IRC_info.get("opt"):
if step_state in (StepState.Reject, StepState.Poor):
step_state = StepState.Reject
if np.isclose(self.trust, params.tmin):
logger.info("IRC stuck with the minimum step-size and bad quality step. Forcing it to take a step.\n")
Expand All @@ -759,7 +764,7 @@ def evaluate_IRC_step(self, params, step_state, criteria_met, IRC_converged):
self.IC_check = True

if self.IRC_info["total_disp"] > 5*self.IRC_init_step:
if criteria_met :
if criteria_met:
if self.IRC_info.get("direction") == 1:
logger.info("\nIRC forward direction converged\n")
logger.info("IRC backward direction starts here\n\n")
Expand All @@ -778,11 +783,17 @@ def evaluate_IRC_step(self, params, step_state, criteria_met, IRC_converged):

return False, step_state

def evaluate_OPT_step(self, params, step_state, criteria_met, Converged_molpro_gmax, Converged_molpro_dmax):
def evaluate_OPT_step(self, params, step_state, criteria_met, Converged_grms, Converged_drms, Converged_energy,
Converged_molpro_gmax, Converged_molpro_dmax):
if criteria_met and self.conSatisfied:
self.SortedEigenvalues(self.H)
logger.info("Converged! =D\n")
self.state = OPT_STATE.CONVERGED
if params.irc and self.IRC_info.get("direction") == 1:
logger.info("\nIRC forward direction converged\n")
logger.info("IRC backward direction starts here\n\n")
self.reset_irc()
else:
self.SortedEigenvalues(self.H)
logger.info("Converged! =D\n")
self.state = OPT_STATE.CONVERGED
return True, step_state

if self.Iteration >= params.maxiter:
Expand Down Expand Up @@ -881,7 +892,8 @@ def evaluateStep(self):
if params.irc and not self.IRC_info.get("opt"):
terminate, step_state = self.evaluate_IRC_step(params, step_state, criteria_met, IRC_converged)
else:
terminate, step_state = self.evaluate_OPT_step(params, step_state, criteria_met, Converged_molpro_gmax, Converged_molpro_dmax)
terminate, step_state = self.evaluate_OPT_step(params, step_state, criteria_met, Converged_grms,
Converged_drms, Converged_energy, Converged_molpro_gmax, Converged_molpro_dmax)

if terminate: return

Expand Down
2 changes: 1 addition & 1 deletion geometric/qcf_neb.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,7 +130,7 @@ def get_basic_info(info_dict, previous=False):
params_dict = {"images": args_dict.get("images"), "maxg": args_dict.get("maximum_force"),
"avgg": args_dict.get("average_force"), "nebk": args_dict.get("spring_constant"),
"neb_maxcyc": args_dict.get("maximum_cycle"), "plain": args_dict.get("spring_type"),
"skip": not args_dict.get("hessian_reset"), "epsilon": args_dict.get("epsilon")}
"epsilon": args_dict.get("epsilon")}
iteration = args_dict.get("iteration")
params = NEBParams(**params_dict)

Expand Down
Loading