from __future__ import print_function
from builtins import super
from datetime import datetime
import numpy as np
from ase.optimize.optimize import Optimizer
from .displacements import get_nvibdof
from .vibrations import harmonic_vibrational_analysis
[docs]class NormalModeBFGS(Optimizer, object):
'''
Normal mode optimizer with approximate hessian update
Parameters
----------
atoms : ase.Atoms
Atoms object with the structure to optimize
phase : str
Phase, should be either `gas` or `solid`
hessian : array_like (N, N)
Initial hessian matrix in eV/Angstrom^2
hessian_update : str
Name of the approximate formula to udpate hessian, one of: `BFGS`,
`SR1`, `DFP`
proj_translations : bool
If ``True`` translational degrees of freedom will be projected from
the hessian
proj_rotations : bool
If ``True`` rotational degrees of freedom will be projected from
the hessian
logfile : str
Name log the log file
trajectory : str
Name of the trajectory file
'''
def __init__(self, atoms, phase, hessian=None,
hessian_update='BFGS',
logfile='-', trajectory=None, restart=None,
proj_translations=True, proj_rotations=True,
master=None, verbose=False):
super().__init__(atoms, restart, logfile, trajectory, master)
self.trajectory = trajectory
self.hessian = hessian
self.phase = phase
self.hessian_update = hessian_update
self.verbose = verbose
self.restart = restart
self.nsteps = 0
self.proj_translations = proj_translations
self.proj_rotations = proj_rotations
# initialize other necessary variables
self.coords_0 = None
self.grad_0 = None
self.natoms = atoms.get_number_of_atoms()
self.ndof = 3 * self.natoms
self.nvibdof = get_nvibdof(atoms, proj_rotations, proj_translations,
self.phase)
# matrix with inverse square roots of masses on diagonal
self.M_invsqrt = np.zeros((self.ndof, self.ndof), dtype=float)
np.fill_diagonal(self.M_invsqrt,
np.repeat(1.0 / np.sqrt(atoms.get_masses()), 3))
# write the header to the logfile
self.log_header()
@property
def hessian(self):
return self._hessian
@hessian.setter
def hessian(self, value):
'Initialize the hessian matrix'
if hessian is None:
self._hessian = np.eye(3 * len(self.atoms)) * 70.0
[docs] def log(self, grad, grad_nm, step_nm):
'Print a line with convergence information'
gmax = np.sqrt((grad.reshape(self.natoms, 3))**2).sum(axis=1).max()
gnnmorm = np.sqrt(np.dot(grad_nm, grad_nm))
print('{0:>6d} {1:>15.8f} {2:>15.8f} {3:>15.8f} {4:>15.8f} {5:>15.8f} {6:>20s}'.format(
self.nsteps, gmax, np.max(np.abs(grad_nm)), gnnmorm,
np.sqrt(np.dot(step_nm, step_nm)), self.atoms.get_potential_energy(),
datetime.now().strftime('%H:%M:%S %d-%m-%Y')), file=self.logfile)
self.logfile.flush()
[docs] def update_hessian(self, coords, grad):
'''
Perform hessian update
Parameters
----------
coords : array_like (N,)
Current coordinates as vector
grad : array_like (N,)
Current gradient
'''
# on first step there's no previous grad and coords so not updated is
# performed
if self.grad_0 is None:
return
macheps = np.finfo(np.float64).eps
dx = coords - self.coords_0
dg = grad - self.grad_0
# same configuration again
if np.abs(dx).max() < 1.0e-7:
return
if self.hessian_update == 'BFGS':
dxdg = np.dot(dx, dg)
hdx = np.dot(self.hessian, dx)
b = np.dot(dx, hdx)
if np.abs(dxdg) > macheps and np.abs(b) > macheps:
self.hessian += np.outer(dg, dg) / dxdg - np.outer(hdx, hdx) / b
elif self.hessian_update == 'DFP':
dgdxT = np.outer(dg, dx)
dgTdx = np.dot(dg, dx)
uleft = np.eye(self.hessian.shape[0]) - dgdxT / dgTdx
uright = np.eye(self.hessian.shape[0]) - dgdxT.T / dgTdx
bk = np.dot(uleft, np.dot(self.hessian, uright))
self.hessian = bk + np.outer(dg, dg) / dgTdx
elif self.hessian_update.upper() == 'SR1':
hdx = np.dot(self.hessian, dx)
dghdx = dg - hdx
self.hessian += np.outer(dghdx, dghdx) / np.dot(dghdx, dx)
else:
raise NotImplementedError('update <{}> not available'.format(self.hessian_update))
[docs] def step(self, grad):
'''
Calculate the step in cartesian coordinates based on the step in
normal modes in the rational function approximation (RFO)
Args:
grad : array_like (N,)
Current gradient
'''
nv = self.nvibdof
coords = self.atoms.get_positions().ravel()
grad = -1.0 * self.atoms.get_forces().ravel()
self.update_hessian(coords, grad)
# calculate hessian eigenvalues and eigenvectors
evals, evecs = harmonic_vibrational_analysis(self.hessian, self.atoms,
proj_translations=self.proj_translations,
proj_rotations=self.proj_rotations,
ascomplex=False, massau=False)
evals = np.power(evals, 2)
mwevecs = np.dot(self.M_invsqrt, evecs)
grad_nm = np.dot(mwevecs.T, grad)
step_nm = np.zeros_like(grad_nm)
step_nm[:nv] = -2.0 * grad_nm[:nv] / (evals[:nv] +
np.sqrt(evals[:nv]**2 + 4.0 * grad_nm[:nv]**2))
step_cart = np.dot(mwevecs, step_nm)
new_coords = coords + step_cart
self.atoms.set_positions(new_coords.reshape(self.natoms, 3))
self.coords_0 = new_coords.copy()
self.grad_0 = grad.copy()
self.dump((self.hessian, self.coords_0, self.grad_0))
self.log(grad, grad_nm, step_nm)
[docs] def run(self, fmax=0.05, steps=100000000):
"""Run structure optimization algorithm.
This method will return when the forces on all individual
atoms are less than *fmax* or when the number of steps exceeds
*steps*."""
self.fmax = fmax
step = 0
while step < steps:
f = self.atoms.get_forces()
self.call_observers()
if self.converged(f):
return
self.step(f)
self.nsteps += 1
step += 1
[docs] def read(self):
self.hessian, self.coords_0, self.grad_0 = self.load()
[docs]def nmoptimize(atoms, hessian, calc, phase, proj_translations=True,
proj_rotations=True, gtol=1.0e-5, verbose=False,
hessian_update='BFGS', steps=100000):
'''
Relax the strcture using normal mode displacements
Parameters
----------
atoms : ase.Atoms
Atoms object with the structure to optimize
hessian : array_like
Hessian matrix in eV/Angstrom^2
calc : ase.Calculator
ASE Calcualtor instance to be used to calculate forces
phase : str
Phase, 'solid' or 'gas'
gtol : float, default=1.0e-5
Energy gradient threshold
hessian_update : str
Approximate formula to update hessian, possible values are 'BFGS',
'SR1' and 'DFP'
steps : int
Maximal number of iteration to be performed
verbose : bool
If ``True`` additional debug information will be printed
Notes
-----
Internally eV and Angstroms are used.
.. seealso::
Bour, P., & Keiderling, T. A. (2002). Partial optimization of molecular geometry
in normal coordinates and use as a tool for simulation of vibrational spectra.
The Journal of Chemical Physics, 117(9), 4126.
`doi:10.1063/1.1498468 <http://dx.doi.org/10.1063/1.1498468>`_
'''
natoms = atoms.get_number_of_atoms()
ndof = 3 * natoms
masses = atoms.get_masses()
coords = atoms.get_positions().ravel()
nvibdof = get_nvibdof(atoms, proj_rotations, proj_translations, phase)
# matrix with inverse square roots of masses on diagonal
M_invsqrt = np.zeros((ndof, ndof), dtype=float)
np.fill_diagonal(M_invsqrt, np.repeat(1.0 / np.sqrt(masses), 3))
# calculate hessian eigenvalues and eigenvectors
evals, evecs = harmonic_vibrational_analysis(hessian, atoms,
proj_translations=proj_translations,
proj_rotations=proj_rotations,
ascomplex=False, massau=False)
evals = np.power(evals, 2)
mwevecs = np.dot(M_invsqrt, evecs)
coords_old = coords.copy()
# run the job for the initial structure
atoms.set_calculator(calc)
# get forces after run
grad = -1.0 * atoms.get_forces().ravel()
grad_old = grad.copy()
grad_nm = np.dot(mwevecs.T, grad)
step_nm = np.zeros_like(grad_nm)
step_nm[:nvibdof] = -2.0 * grad_nm[:nvibdof] / (evals[:nvibdof]
+ np.sqrt(evals[:nvibdof]**2 + 4.0 * grad_nm[:nvibdof]**2))
step_cart = np.dot(mwevecs, step_nm)
coords = coords_old + step_cart
if verbose:
print(' eigenvalues '.center(50, '-'))
print(evals)
print(' cart gradient '.center(50, '-'))
print(grad)
print(' nm gradient '.center(50, '-'))
print(grad_nm)
print(' nm step '.center(50, '-'))
print(step_nm)
print(' cart step '.center(50, '-'))
print(step_cart)
print(' new coordinates '.center(50, '-'))
print(coords)
iteration = 0
# header for the convergence information
print('{0:<6s} {1:^15s} {2:^15s} {3:^15s} {4:^15s} {5:^20s}'.format(
'iter', 'G(NM) max', 'G(NM) norm', 'NM step norm', 'energy [eV]', 'time'))
print('=' * 91)
print('{0:>6d} {1:>15.8f} {2:>15.8f} {3:>15.8f} {4:>15.8f} {5:>20s}'.format(
iteration, np.max(np.abs(grad_nm)), np.sqrt(np.dot(grad_nm, grad_nm)),
np.sqrt(np.dot(step_nm, step_nm)), atoms.get_potential_energy(),
datetime.now().strftime('%H:%M:%S %d-%m-%Y')))
while iteration <= steps:
iteration += 1
# delta_coord = coords - coords_old
atoms.set_positions(coords.reshape(natoms, 3))
coords_old = coords.copy()
grad = -1.0 * atoms.get_forces().ravel()
# delta_grad = grad - grad_old
hessian = update_hessian(grad, grad_old, step_cart, hessian,
update=hessian_update)
grad_old = grad.copy()
# calculate hessian eigenvalues and eigenvectors
evals, evecs = harmonic_vibrational_analysis(hessian, atoms,
proj_translations=proj_translations,
proj_rotations=proj_rotations,
ascomplex=False, massau=False)
evals = np.power(evals, 2)
mwevecs = np.dot(M_invsqrt, evecs)
grad_nm = np.dot(mwevecs.T, grad)
gmax = np.max(np.abs(grad_nm))
# print the convergence info
print('{0:>6d} {1:>15.8f} {2:>15.8f} {3:>15.8f} {4:>15.8f} {5:>20s}'.format(
iteration, gmax, np.sqrt(np.dot(grad_nm, grad_nm)),
np.sqrt(np.dot(step_nm, step_nm)), atoms.get_potential_energy(),
datetime.now().strftime('%H:%M:%S %d-%m-%Y')))
if gmax < gtol:
evals, evecs = harmonic_vibrational_analysis(hessian, atoms,
proj_translations=proj_translations,
proj_rotations=proj_rotations,
ascomplex=False, massau=False)
np.save('hessian_evalues', evals)
np.save('hessian_evectors', evecs)
print('# convergence achieved after {} iterations'.format(iteration))
break
step_nm[:nvibdof] = -2.0 * grad_nm[:nvibdof] / (evals[:nvibdof]
+ np.sqrt(evals[:nvibdof]**2 + 4.0 * grad_nm[:nvibdof]**2))
step_cart = np.dot(mwevecs, step_nm)
coords = coords_old + step_cart
if verbose:
print(' eigenvalues '.center(50, '-'))
print(evals)
print(' cart gradient '.center(50, '-'))
print(grad)
print(' nm gradient '.center(50, '-'))
print(grad_nm)
print(' nm step '.center(50, '-'))
print(step_nm)
print(' cart step '.center(50, '-'))
print(step_cart)
print(' new coordinates '.center(50, '-'))
print(coords)
else:
print('### convergence NOT achieved after ', iteration, ' iterations')
[docs]def update_hessian(grad, grad_old, dx, hessian, update='BFGS'):
'''
Perform hessian update
Parameters
----------
grad : array_like (N,)
Current gradient
grad_old : array_like (N,)
Previous gradient
dx : array_like (N,)
Step vector x_n - x_(n-1)
hessian : array_like (N, N)
Hessian matrix
update : str
Name of the hessian update to perform, possible values are 'BFGS',
'SR1' and 'DFP'
Returns
-------
hessian : array_like
Update hessian matrix
'''
macheps = np.finfo(np.float64).eps
dg = grad - grad_old
if update == 'BFGS':
dxdg = np.dot(dx, dg)
hdx = np.dot(hessian, dx)
b = np.dot(dx, hdx)
if np.abs(dxdg) < macheps or np.abs(b) < macheps:
return hessian
else:
return hessian + np.outer(dg, dg) / dxdg - np.outer(hdx, hdx) / b
elif update == 'DFP':
dgdxT = np.outer(dg, dx)
dgTdx = np.dot(dg, dx)
uleft = np.eye(hessian.shape[0]) - dgdxT / dgTdx
uright = np.eye(hessian.shape[0]) - dgdxT.T / dgTdx
bk = np.dot(uleft, np.dot(hessian, uright))
return bk + np.outer(dg, dg) / dgTdx
elif update.upper() == 'SR1':
hdx = np.dot(hessian, dx)
dghdx = dg - hdx
return hessian + np.outer(dghdx, dghdx) / np.dot(dghdx, dx)
else:
raise NotImplementedError