Skip to content

Commit

Permalink
Merge pull request #207 from hjnpark/irc_update
Browse files Browse the repository at this point in the history 
IRC bug fix and NEB update
  • Loading branch information
@leeping
leeping authored Jul 19, 2024
2 parents 740fcec + d92ef1b commit 00e7722
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 17 deletions.
44 changes: 28 additions & 16 deletions geometric/optimize.py
View file
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
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

0 comments on commit 00e7722

Please sign in to comment.