Commit 845a2899 authored by Christopher Ostrouchov's avatar Christopher Ostrouchov
Browse files

Basic implementation of dftfit

parent b4052730
from itertols import combinations
from itertools import combinations
from tempfile import TemporaryDirectory
import math
import numpy as np
from .io.lammps import LammpsWriter, LammpsRunner
......@@ -22,7 +25,7 @@ def optimize_function(md_calculations, dft_calculations, weight_forces, weight_s
n_energy_sq_error = 0.0
d_energy_sq_error = 0.0
for md_calculation, dft_calcualtion in zip(md_calculations, dft_calculations):
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)
......@@ -33,11 +36,11 @@ def optimize_function(md_calculations, dft_calculations, weight_forces, weight_s
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 = sqrt(n_force_sq_error / d_force_sq_error)
stress_sq_error = sqrt(n_stress_sq_error / d_stress_sq_error)
energy_sq_error = sqrt(n_energy_sq_error / d_energy_sq_error)
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)
return (weight_force * force_sq_error +
return (weight_forces * force_sq_error +
weight_stress * stress_sq_error +
weight_energy * energy_sq_error)
......
import numpy as np
from scipy import optimize
from ._optimize import evaluate, optimize_function
class Dftfit:
......@@ -11,26 +14,55 @@ class Dftfit:
MD_SOLVERS = {'LAMMPS'}
def __init__(self, max_iter=-1, step_toll=1e-6, w_f=0.8, w_s=0.1, w_e=0.1):
self.max_iter = max_iter
self.step_toll = step_toll
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"):
self.max_iterations = max_iterations
self.step_tolerance = step_tolerance
self.method = method
if not np.isclose(sum([w_f, w_s, w_e]), 1.0, 1e-8):
if not np.isclose(sum([weight_forces, weight_stress, weight_energy]), 1.0, 1e-8):
raise ValueError('sum of weighting functions must be one')
self.weights = {
'forces': w_f,
'stress': w_s,
'energy': w_e
'forces': weight_forces,
'stress': weight_stress,
'energy': weight_energy
}
def fit(self, dft_calculations, starting_potential):
def fit(self, calculations, potential):
"""
Args:
- calculations: training set of DFT calculations to fit
- potential: model of Potential to fit
"""
print('Fitting.....')
def optimization_function(parameters):
potential.optimization_parameters = parameters
md_calculations = []
for calculation in calculations:
md_calculation = evaluate('LAMMPS', calculation.structure, potential)
md_calculations.append(md_calculation)
score = optimize_function(md_calculations, calculations,
weight_forces=self.weights['forces'],
weight_stress=self.weights['stress'],
weight_energy=self.weights['energy'])
print(' '.join(['{:12.8g}'.format(_) for _ in parameters]), score)
return score
result = optimize.minimize(
fun=optimization_function,
x0=potential.optimization_parameters,
method=self.method,
bounds=potential.optimization_bounds,
tol=self.step_tolerance,
options={'disp': True, 'maxiter': self.max_iterations}
)
print(result)
print('Done!')
def predict(self, structure):
pass
......@@ -5,16 +5,21 @@ class FloatParameter:
""" Float with tracking. initial value and bounds.
"""
def __init__(self, initial, bounds=(-sys.float_info.max, sys.float_info.max), computed=None):
self.initial = float(initial)
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.computed = computed
self.computed = None
self._fixed = fixed
@property
def fixed(self):
return self._fixed or self.computed != None
def __float__(self):
if self.computed is None:
return self.current
return self.computed()
if self.computed:
return self.computed()
return self.current
def __str__(self):
return str(self.current)
......@@ -27,7 +27,7 @@ class Potential:
if not {e.symbol for e in composition.keys()} <= charges.keys():
raise ValueError('charge ballance constrains requires all elements to be defined in charge')
for charge_element, parameter in charges.items():
if isinstance(parameter, FloatParameter) and parameter.computed == None:
if isinstance(parameter, FloatParameter) and not parameter.fixed:
break
else:
if abs(sum(float(charges[element.symbol]) * amount for element, amount in composition.items())) > 1e-8:
......@@ -38,17 +38,19 @@ class Potential:
def _collect_parameters(self):
self._parameters = []
def _walk(value):
def _walk(value): # Ordered traversal of dictionary
if isinstance(value, dict):
for key in value:
for key in sorted(value.keys()):
_walk(value[key])
elif isinstance(value, (tuple, list)):
for item in value:
_walk(item)
elif isinstance(value, FloatParameter) and value.computed == None:
elif isinstance(value, FloatParameter):
self._parameters.append(value)
_walk(self.schema)
self._optimization_parameters = [p for p in self._parameters if not p.fixed]
@classmethod
def from_file(cls, filename, format=None):
if format not in {'json', 'yaml'}:
......@@ -73,16 +75,24 @@ class Potential:
"""
return np.array([float(parameter) for parameter in self._parameters])
@parameters.setter
def parameters(self, parameters):
@property
def optimization_parameters(self):
return self._optimization_parameters
@optimization_parameters.setter
def optimization_parameters(self, parameters):
""" Update potential with given parameters
"""
if len(parameters) != len(self._parameters):
if len(parameters) != len(self._optimization_parameters):
raise ValueError('updating parameters does not match length of potential parameters')
for parameter, update_parameter in zip(self._parameters, parameters):
for parameter, update_parameter in zip(self._optimization_parameters, parameters):
parameter.current = update_parameter
@property
def optimization_bounds(self):
return np.array([parameter.bounds for parameter in self._optimization_parameters])
def __str__(self):
return json.dumps(self.schema, sort_keys=True, indent=4)
......@@ -12,10 +12,11 @@ class ParameterSchema(BaseSchema):
bounds = fields.List(fields.Float(), validate=validate.Length(equal=2), missing=(-sys.float_info.max, sys.float_info.max))
class FloatOrParameter(fields.Field):
class FloatParameterField(fields.Field):
def _deserialize(self, value, attr, data):
try:
return float(value)
value = float(value)
return FloatParameter(value, fixed=True)
except (ValueError, TypeError):
schema_load, errors = ParameterSchema().load(value)
return FloatParameter(**value)
......@@ -36,26 +37,26 @@ class ConstraintSchema(BaseSchema):
charge_balance = fields.String()
ChargesSchema = type('ChargeSchema', (BaseSchema,), {element: FloatOrParameter() for element in element_symbols})
ChargesSchema = type('ChargeSchema', (BaseSchema,), {element: FloatParameterField() for element in element_symbols})
class KspaceSchema(BaseSchema):
KSPACE_TYPES = {'ewald', 'pppm'}
type = fields.String(required=True, validate=validate.OneOf(KSPACE_TYPES))
tollerance = FloatOrParameter(required=True, validate=validate.Range(min=1e-16))
tollerance = FloatParameterField(required=True, validate=validate.Range(min=1e-16))
class ParametersSchema(BaseSchema):
elements = fields.List(fields.String(validate=validate.OneOf(element_symbols)), required=True)
coefficients = fields.List(FloatOrParameter())
coefficients = fields.List(FloatParameterField())
class PairPotentialSchema(BaseSchema):
PAIR_POTENTIALS = {'buckingham'}
type = fields.String(required=True, validate=validate.OneOf(PAIR_POTENTIALS))
cutoff = FloatOrParameter(required=True, validate=validate.Range(min=1e-6))
cutoff = FloatParameterField(required=True, validate=validate.Range(min=1e-6))
parameters = fields.Nested(ParametersSchema, required=True, many=True)
......
......@@ -26,6 +26,9 @@ class Training:
def calculations(self):
return self._calculations
def __iter__(self):
return iter(self._calculations)
def download_mattoolkit_calculations(self, selector):
from mattoolkit.api import CalculationResourceList
......
......@@ -33,7 +33,7 @@ with open(path.join(here, 'README.md'), encoding='utf-8') as f:
long_description = f.read()
version='0.0.1'
version='0.0.2'
setup(
name='dftfit',
version=version,
......
......@@ -32,4 +32,5 @@ def test_potential_from_file():
def test_potential_from_file_read_parameters():
potential = Potential.from_file('test_files/potential/mgo-fitting.yaml')
assert len(potential.parameters) == 8
assert len(potential.parameters) == 13
assert len(potential.optimization_parameters) == 8
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