Commit 2ef9d0e4 authored by Christopher Ostrouchov's avatar Christopher Ostrouchov
Browse files

Normalization of parameters

parent 751dc814
from itertools import combinations
from tempfile import TemporaryDirectory
import math
import numpy as np
from .io.lammps import LammpsWriter, LammpsRunner
def evaluate(runner, structure, potential):
if runner == 'LAMMPS':
with TemporaryDirectory() as tempdir:
writer = LammpsWriter(structure, potential)
return LammpsRunner().run(writer, ['lammps', '-i', 'lammps.in'], tempdir)
else:
raise ValueError('unknown runner')
def optimize_function(md_calculations, dft_calculations, weights):
n_force_sq_error = 0.0
d_force_sq_error = 0.0
n_stress_sq_error = 0.0
d_stress_sq_error = 0.0
n_energy_sq_error = 0.0
d_energy_sq_error = 0.0
for md_calculation, dft_calculation in zip(md_calculations, dft_calculations):
n_force_sq_error += np.sum((md_calculation.forces - dft_calculation.forces)**2.0)
d_force_sq_error += np.sum(dft_calculation.forces**2.0)
n_stress_sq_error += np.sum((md_calculation.stress - dft_calculation.stress)**2.0)
d_stress_sq_error += np.sum(dft_calculation.stress**2.0)
for (md_calc_i, dft_calc_i), (md_calc_j, dft_calc_j) in combinations(zip(md_calculations, dft_calculations), 2):
n_energy_sq_error += ((md_calc_i.energy - md_calc_j.energy) - (dft_calc_i.energy - dft_calc_j.energy))**2.0
d_energy_sq_error += (dft_calc_i.energy - dft_calc_j.energy)**2.0
force_sq_error = math.sqrt(n_force_sq_error / d_force_sq_error)
stress_sq_error = math.sqrt(n_stress_sq_error / d_stress_sq_error)
energy_sq_error = math.sqrt(n_energy_sq_error / d_energy_sq_error)
score = (
weights['forces'] * force_sq_error + \
weights['stress'] * stress_sq_error + \
weights['energy'] * energy_sq_error
)
return {
'weights': {'forces': weights['forces'], 'stress': weights['stress'], 'energy': weights['energy']},
'parts': {'forces': force_sq_error, 'stress': stress_sq_error, 'energy': energy_sq_error},
'score': score
}
# def optimize():
# # where to optimize potential based on condition
# if package == 'scipy':
# minimize(
# fun=optimize_function,
# x0=self.initial_parameters,
# method="L-BFGS-B",
# #method="powell",
# #method="COBYLA",
# #method="TNC",
# #method="SLSQP",
# jac=False,
# bounds=self.bounds,
# options={'disp': True},
# )
# elif package == 'nlopt':
# # NLOPT optimization solver
# # used because reference exists that shows algorithm works
# # Using BOBYQA algorithm
# # http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms#BOBYQA
# opt = nlopt.opt(nlopt.LN_BOBYQA, len(self.initial_parameters))
# lower_bounds, upper_bounds = zip(*self.bounds)
# self.logger.info((
# "Using {} optimization algorithm with {} parameters\n"
# "{}\n"
# ).format(opt.get_algorithm(), opt.get_dimension(), opt.get_algorithm_name()))
# opt.set_lower_bounds(lower_bounds)
# opt.set_upper_bounds(upper_bounds)
# self.logger.info((
# "\nBounded Problem:\n"
# "Lower Bound: {}\n"
# "Upper Bound: {}\n"
# ).format(opt.get_lower_bounds(), opt.get_upper_bounds()))
# opt.set_min_objective(optimize_function)
# opt.set_ftol_rel(1e-6)
# # opt.set_maxeval(0) # 0 or negative for no limit
# # opt.set_maxtime(0) # seconds (0 or negative for no limit)
# optimized_parameters = opt.optimize(list(self.initial_parameters))
# min_objective_value = opt.last_optimum_value()
# else:
# raise ValueError('Optimzation Packages %s not recognized' % package)
# pass
import numpy as np
from scipy import optimize
from ._optimize import evaluate, optimize_function
from .optimize import evaluate, optimize_function
class Dftfit:
......@@ -16,9 +16,11 @@ class Dftfit:
def __init__(self,
weight_forces=0.8, weight_stress=0.1, weight_energy=0.1,
max_iterations=100, step_tolerance=1e-6, method="L-BFGS-B"):
step_jacobian=1e-5, method="L-BFGS-B",
max_iterations=100, step_tolerance=1e-6):
self.max_iterations = max_iterations
self.step_tolerance = step_tolerance
self.step_jacobian = step_jacobian
self.method = method
if not np.isclose(sum([weight_forces, weight_stress, weight_energy]), 1.0, 1e-8):
......@@ -30,6 +32,15 @@ class Dftfit:
'energy': weight_energy
}
def _normalize_parameters(self, potential):
bound_range = potential.optimization_bounds[:, 1] - potential.optimization_bounds[:, 0]
range_within_limit = (1e-8 < bound_range) & (bound_range < 1e12)
subtraction_vector = np.where(range_within_limit, potential.optimization_bounds[:, 0], 0)
scale_vector = np.where(range_within_limit, potential.optimization_bounds[:, 1], 1)
self._normalized_bounds = np.where(np.column_stack((range_within_limit, range_within_limit)), np.array([0, 1]), potential.optimization_bounds)
self._normalize_vector = lambda v: (v - subtraction_vector) / scale_vector
self._unnormalize_vector = lambda v: (v * scale_vector) + subtraction_vector
def fit(self, calculations, initial_potential):
"""
......@@ -38,9 +49,12 @@ class Dftfit:
- calculations: training set of DFT calculations to fit
- potential: model of Potential to fit
"""
self._normalize_parameters(initial_potential)
initial_normalized_parameters = self._normalize_vector(initial_potential.optimization_parameters)
def optimization_function(parameters):
potential = initial_potential.copy()
potential.optimization_parameters = parameters
potential.optimization_parameters = self._unnormalize_vector(parameters)
md_calculations = []
for calculation in calculations:
......@@ -48,18 +62,24 @@ class Dftfit:
md_calculations.append(md_calculation)
result = optimize_function(md_calculations, calculations, self.weights)
parameter_str = ' '.join(['{:12.8g}'.format(_) for _ in parameters])
parameter_str = ' '.join(['{:12.8g}'.format(_) for _ in (potential.optimization_parameters)])
print(parameter_str)
parameter_str = ' '.join(['{:12.8g}'.format(_) for _ in (potential.optimization_parameters - initial_potential.optimization_parameters)])
optimization_str = '%12.6f %12.6f %12.6f %12.6f' % (result['parts']['forces'], result['parts']['stress'], result['parts']['energy'], result['score'])
print(parameter_str, '|', optimization_str)
return result['score']
print('It is recommended to bound all parameters for better convergence')
print('Initial Parameters:')
print(' '.join(['{:12.8g}'.format(_) for _ in initial_potential.optimization_parameters]))
result = optimize.minimize(
fun=optimization_function,
x0=initial_potential.optimization_parameters,
x0=initial_normalized_parameters,
method=self.method,
bounds=initial_potential.optimization_bounds,
bounds=self._normalized_bounds,
jac=False,
tol=self.step_tolerance,
options={'disp': True, 'maxiter': self.max_iterations}
options={'disp': True, 'maxiter': self.max_iterations, 'eps': self.step_jacobian}
)
print(result)
......
This diff is collapsed.
......@@ -8,7 +8,7 @@ class FloatParameter:
def __init__(self, initial, bounds=None, fixed=False):
bounds = bounds or (-sys.float_info.max, sys.float_info.max)
self.current = float(initial)
self.bounds = [float(_) for _ in bounds]
self._bounds = [float(_) for _ in bounds]
self.computed = None
self._fixed = fixed
......@@ -16,6 +16,10 @@ class FloatParameter:
def fixed(self):
return self._fixed or self.computed != None
@property
def bounds(self):
return self._bounds
def __float__(self):
if self.computed:
return self.computed()
......
......@@ -49,7 +49,12 @@ class Potential:
self._parameters.append(value)
_walk(self.schema)
self._optimization_parameters = [p for p in self._parameters if not p.fixed]
self._optimization_parameters = []
self._optimization_parameter_indicies = []
for i, p in enumerate(self._parameters):
if not p.fixed:
self._optimization_parameters.append(p)
self._optimization_parameter_indicies.append(i)
@classmethod
def from_file(cls, filename, format=None):
......@@ -116,6 +121,10 @@ class Potential:
def optimization_bounds(self):
return np.array([parameter.bounds for parameter in self._optimization_parameters])
@property
def optimization_parameter_indicies(self):
return self._optimization_parameter_indicies
def __repr__(self):
return self.__str__()
......
......@@ -50,7 +50,7 @@ setup(
keywords='materials dft molecular dynamics lammps science hpc',
download_url='https://gitlab.aves.io/costrouc/dftfit/repository/archive.zip?ref=v%s' % version,
packages=find_packages(exclude=('tests', 'docs', 'notebooks', 'examples')),
install_requires=['pymatgen==2017.7.4', 'pymatgen-lammps', 'marshmallow', 'pyyaml', 'click'],
install_requires=['pymatgen==2017.7.4', 'pymatgen-lammps', 'marshmallow', 'pyyaml', 'click', 'sqlalchemy'],
extras_require={
'mattoolkit': 'mattoolkit'
},
......
Supports Markdown
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