PK!#  berny/Logger.py# This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. from __future__ import print_function import re class Logger(object): """ Manages logging of information. :param int verbosity: verbosity level :param file out: argument passed to ``print(..., file)`` :param str regex: log only messages that match regex """ def __init__(self, verbosity=None, out=None, regex=None): self.verbosity = verbosity or 0 self.out = out self.n = 0 self.regex = regex def __call__(self, msg, level=0): """ Log a message. :param str msg: the message :param int level: this is compared against the verbosity level """ if level < -self.verbosity: return if self.regex is not None and not re.search(self.regex, msg): return print(self.n, msg, file=self.out) PK!w( berny/Math.py# This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. import numpy as np def rms(A): if A.size == 0: return None return np.sqrt(np.sum(A**2)/A.size) def pinv(A, log=lambda _: None): U, D, V = np.linalg.svd(A) thre = 1e3 thre_log = 1e8 gaps = D[:-1]/D[1:] try: n = np.flatnonzero(gaps > thre)[0] except IndexError: n = len(gaps) else: gap = gaps[n] if gap < thre_log: log('Pseudoinverse gap of only: {:.1e}'.format(gap)) D[n+1:] = 0 D[:n+1] = 1/D[:n+1] return U.dot(np.diag(D)).dot(V) def cross(a, b): return np.array([ a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0] ]) def fit_cubic(y0, y1, g0, g1): """Fit cubic polynomial to function values and derivatives at x = 0, 1. Returns position and function value of minimum if fit succeeds. Fit does not succeeds if 1. polynomial doesn't have extrema or 2. maximum is from (0,1) or 3. maximum is closer to 0.5 than minimum """ a = 2*(y0-y1)+g0+g1 b = -3*(y0-y1)-2*g0-g1 p = np.array([a, b, g0, y0]) r = np.roots(np.polyder(p)) if not np.isreal(r).all(): return None, None r = sorted(x.real for x in r) if p[0] > 0: maxim, minim = r else: minim, maxim = r if 0 < maxim < 1 and abs(minim-0.5) > abs(maxim-0.5): return None, None return minim, np.polyval(p, minim) def fit_quartic(y0, y1, g0, g1): """Fit constrained quartic polynomial to function values and erivatives at x = 0,1. Returns position and function value of minimum or None if fit fails or has a maximum. Quartic polynomial is constrained such that it's 2nd derivative is zero at just one point. This ensures that it has just one local extremum. No such or two such quartic polynomials always exist. From the two, the one with lower minimum is chosen. """ def g(y0, y1, g0, g1, c): a = c+3*(y0-y1)+2*g0+g1 b = -2*c-4*(y0-y1)-3*g0-g1 return np.array([a, b, c, g0, y0]) def quart_min(p): r = np.roots(np.polyder(p)) is_real = np.isreal(r) if is_real.sum() == 1: minim = r[is_real][0].real else: minim = r[(r == max(-abs(r))) | (r == -max(-abs(r)))][0].real return minim, np.polyval(p, minim) D = -(g0+g1)**2-2*g0*g1+6*(y1-y0)*(g0+g1)-6*(y1-y0)**2 # discriminant of d^2y/dx^2=0 if D < 1e-11: return None, None else: m = -5*g0-g1-6*y0+6*y1 p1 = g(y0, y1, g0, g1, .5*(m+np.sqrt(2*D))) p2 = g(y0, y1, g0, g1, .5*(m-np.sqrt(2*D))) if p1[0] < 0 and p2[0] < 0: return None, None [minim1, minval1] = quart_min(p1) [minim2, minval2] = quart_min(p2) if minval1 < minval2: return minim1, minval1 else: return minim2, minval2 class FindrootException(Exception): pass def findroot(f, lim): """Find root of increasing function on (-inf,lim). Assumes f(-inf) < 0, f(lim) > 0. """ d = 1. for _ in range(1000): val = f(lim-d) if val > 0: break d = d/2 # find d so that f(lim-d) > 0 else: raise RuntimeError('Cannot find f(x) > 0') x = lim-d # initial guess dx = 1e-10 # step for numerical derivative fx = f(x) err = abs(fx) for _ in range(1000): fxpdx = f(x+dx) dxf = (fxpdx-fx)/dx x = x-fx/dxf fx = f(x) err_new = abs(fx) if err_new >= err: return x err = err_new else: raise FindrootException() PK!&pberny/__init__.pyfrom .optimizers import optimize from .berny import Berny from .Logger import Logger from . import geomlib from .geomlib import Geometry from .coords import angstrom PK!vYqXX berny/_py2.py# taken from Python 3 stdlib class Generator(object): def __iter__(self): return self def next(self): return self.__next__() def throw(self, typ, val=None, tb=None): if val is None: if tb is None: raise typ val = typ() if tb is not None: val = val.with_traceback(tb) raise val def close(self): try: self.throw(GeneratorExit) except (GeneratorExit, StopIteration): pass else: raise RuntimeError("generator ignored GeneratorExit") PK!T2%%berny/berny.py# This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. import sys from collections import namedtuple from itertools import chain import numpy as np from numpy import dot, eye from numpy.linalg import norm from . import Math from .coords import InternalCoords from .Logger import Logger if sys.version_info[:2] >= (3, 5): from collections.abc import Generator else: from ._py2 import Generator # noqa __version__ = '0.2.1' defaults = { 'gradientmax': 0.45e-3, 'gradientrms': 0.3e-3, 'stepmax': 1.8e-3, 'steprms': 1.2e-3, 'trust': 0.3, 'dihedral': True, 'superweakdih': False, } """ - gradientmax, gradientrms, stepmax, steprms: Convergence criteria in atomic units ("step" refers to the step in internal coordinates, assuming radian units for angles). - trust: Initial trust radius in atomic units. It is the maximum RMS of the quadratic step (see below). - dihedral: Form dihedral angles. - superweakdih: Form dihedral angles containing two or more noncovalent bonds. """ class Berny(Generator): """ Generator that receives energy and gradients and yields the next geometry. :param Gometry geom: geometry to start with :param Logger log: used for logging if given :param bool debug: if True, the generator yields debug info on receiving the energy and gradients, otherwise it yields None :param dict restart: start from a state saved from previous run using ``debug=True`` :param int maxsteps: abort after maximum number of steps :param int verbosity: if present and log is None, specifies the verbosity of the default :py:class:`~berny.Logger` :param params: parameters that override the :py:data:`~berny.berny.defaults` The Berny object is to be used as follows:: optimizer = Berny(geom) for geom in optimizer: # calculate energy and gradients (as N-by-3 matrix) debug = optimizer.send((energy, gradients)) """ class State(object): pass Point = namedtuple('Point', 'q E g') def __init__(self, geom, log=None, debug=False, restart=None, maxsteps=100, verbosity=None, **params): self._log = log or Logger(verbosity=verbosity or 0) self._debug = debug self._maxsteps = maxsteps self._converged = False self._n = 0 s = self._state = Berny.State() if restart: vars(s).update(restart) return s.geom = geom s.params = dict(chain(defaults.items(), params.items())) s.trust = s.params['trust'] s.coords = InternalCoords( s.geom, dihedral=s.params['dihedral'], superweakdih=s.params['superweakdih'], ) s.H = s.coords.hessian_guess(s.geom) s.weights = s.coords.weights(s.geom) for line in str(s.coords).split('\n'): self._log(line) s.future = Berny.Point(s.coords.eval_geom(s.geom), None, None) s.first = True def __next__(self): assert self._n <= self._maxsteps if self._n == self._maxsteps or self._converged: raise StopIteration self._n += 1 return self._state.geom def send(self, energy_gradients): log = self._log log.n = self._n s = self._state energy, gradients = energy_gradients gradients = np.array(gradients) log('Energy: {:.12}'.format(energy), level=1) B = s.coords.B_matrix(s.geom) B_inv = B.T.dot(Math.pinv(np.dot(B, B.T), log=log)) current = Berny.Point( s.future.q, energy, dot(B_inv.T, gradients.reshape(-1)) ) if not s.first: s.H = update_hessian( s.H, current.q-s.best.q, current.g-s.best.g, log=log ) s.trust = update_trust( s.trust, current.E-s.previous.E, s.predicted.E-s.interpolated.E, s.predicted.q-s.interpolated.q, log=log ) dq = s.best.q-current.q t, E = linear_search( current.E, s.best.E, dot(current.g, dq), dot(s.best.g, dq), log=log ) s.interpolated = Berny.Point( current.q+t*dq, E, t*s.best.g+(1-t)*current.g ) else: s.interpolated = current if s.trust < 1e-6: raise RuntimeError('The trust radius got too small, check forces?') proj = dot(B, B_inv) H_proj = proj.dot(s.H).dot(proj) + 1000*(eye(len(s.coords))-proj) dq, dE, on_sphere = quadratic_step( dot(proj, s.interpolated.g), H_proj, s.weights, s.trust, log=log ) s.predicted = Berny.Point(s.interpolated.q+dq, s.interpolated.E+dE, None) dq = s.predicted.q-current.q log('Total step: RMS: {:.3}, max: {:.3}'.format( Math.rms(dq), max(abs(dq)) )) q, s.geom = s.coords.update_geom( s.geom, current.q, s.predicted.q-current.q, B_inv, log=log ) s.future = Berny.Point(q, None, None) s.previous = current if s.first or current.E < s.best.E: s.best = current s.first = False self._converged = is_converged( gradients, s.future.q-current.q, on_sphere, s.params, log=log ) if self._n == self._maxsteps: log('Maximum number of steps reached') if self._debug: return vars(s).copy() def throw(self, *args, **kwargs): return Generator.close(self, *args, **kwargs) def no_log(msg, **kwargs): pass def update_hessian(H, dq, dg, log=no_log): dH = dg[None, :]*dg[:, None]/dot(dq, dg) - \ H.dot(dq[None, :]*dq[:, None]).dot(H)/dq.dot(H).dot(dq) # BFGS update log('Hessian update information:') log('* Change: RMS: {:.3}, max: {:.3}'.format(Math.rms(dH), abs(dH).max())) return H+dH def update_trust(trust, dE, dE_predicted, dq, log=no_log): if dE != 0: r = dE/dE_predicted # Fletcher's parameter else: r = 1. log("Trust update: Fletcher's parameter: {:.3}".format(r)) if r < 0.25: return norm(dq)/4 elif r > 0.75 and abs(norm(dq)-trust) < 1e-10: return 2*trust else: return trust def linear_search(E0, E1, g0, g1, log=no_log): log('Linear interpolation:') log('* Energies: {:.8}, {:.8}'.format(E0, E1)) log('* Derivatives: {:.3}, {:.3}'.format(g0, g1)) t, E = Math.fit_quartic(E0, E1, g0, g1) if t is None or t < -1 or t > 2: t, E = Math.fit_cubic(E0, E1, g0, g1) if t is None or t < 0 or t > 1: if E0 <= E1: log('* No fit succeeded, staying in new point') return 0, E0 else: log('* No fit succeeded, returning to best point') return 1, E1 else: msg = 'Cubic interpolation was performed' else: msg = 'Quartic interpolation was performed' log('* {}: t = {:.3}'.format(msg, t)) log('* Interpolated energy: {:.8}'.format(E)) return t, E def quadratic_step(g, H, w, trust, log=no_log): ev = np.linalg.eigvalsh((H+H.T)/2) rfo = np.vstack((np.hstack((H, g[:, None])), np.hstack((g, 0))[None, :])) D, V = np.linalg.eigh((rfo+rfo.T)/2) dq = V[:-1, 0]/V[-1, 0] l = D[0] if norm(dq) <= trust: log('Pure RFO step was performed:') on_sphere = False else: def steplength(l): return norm(np.linalg.solve(l*eye(H.shape[0])-H, g))-trust l = Math.findroot(steplength, ev[0]) # minimization on sphere dq = np.linalg.solve(l*eye(H.shape[0])-H, g) on_sphere = True log('Minimization on sphere was performed:') dE = dot(g, dq)+0.5*dq.dot(H).dot(dq) # predicted energy change log('* Trust radius: {:.2}'.format(trust)) log('* Number of negative eigenvalues: {}'.format((ev < 0).sum())) log('* Lowest eigenvalue: {:.3}'.format(ev[0])) log('* lambda: {:.3}'.format(l)) log('Quadratic step: RMS: {:.3}, max: {:.3}'.format(Math.rms(dq), max(abs(dq)))) log('* Predicted energy change: {:.3}'.format(dE)) return dq, dE, on_sphere def is_converged(forces, step, on_sphere, params, log=no_log): criteria = [ ('Gradient RMS', Math.rms(forces), params['gradientrms']), ('Gradient maximum', np.max(abs(forces)), params['gradientmax']) ] if on_sphere: criteria.append(('Minimization on sphere', False)) else: criteria.extend([ ('Step RMS', Math.rms(step), params['steprms']), ('Step maximum', np.max(abs(step)), params['stepmax']) ]) log('Convergence criteria:') all_matched = True for crit in criteria: if len(crit) > 2: result = crit[1] < crit[2] msg = '{:.3} {} {:.3}'.format(crit[1], '<' if result else '>', crit[2]) else: msg, result = crit msg = '{}: {}'.format(crit[0], msg) if msg else crit[0] msg = '* {} => {}'.format(msg, 'OK' if result else 'no') log(msg) if not result: all_matched = False if all_matched: log('* All criteria matched', level=1) return all_matched PK!JPB B berny/cli.py# This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. import sys import json import pickle from argparse import ArgumentParser from socket import socket, AF_INET, SOCK_STREAM from contextlib import contextmanager from berny import Berny, geomlib @contextmanager def berny_unpickled(berny=None): picklefile = 'berny.pickle' if not berny: with open(picklefile, 'rb') as f: berny = pickle.load(f) try: yield berny except: raise with open(picklefile, 'wb') as f: pickle.dump(berny, f) def handler(berny, f): energy = float(next(f)) gradients = [[float(x) for x in l.split()] for l in f if l.strip()] berny.send(energy, gradients) return next(berny) def get_berny(args): geom = geomlib.load(sys.stdin, args.format) if args.paramfile: with open(args.paramfile) as f: params = json.load(f) else: params = None berny = Berny(geom, **params) next(berny) return berny def init(args): berny = get_berny(args) berny.geom_format = args.format with berny_unpickled(berny) as berny: pass def server(args): berny = get_berny(args) host, port = args.socket server = socket(AF_INET, SOCK_STREAM) server.bind((host, int(port))) server.listen(0) while True: sock, addr = server.accept() f = sock.makefile('r+') geom = handler(berny, f) if geom: f.write(geom.dumps(args.format)) f.flush() f.close() sock.close() if not geom: break def driver(): try: with berny_unpickled() as berny: geom = handler(berny, sys.stdin) except FileNotFoundError: sys.stderr.write('error: No pickled berny, run with --init first?\n') sys.exit(1) if not geom: sys.exit(0) geom.dump(sys.stdout, berny.geom_format) sys.exit(10) def main(): parser = ArgumentParser() arg = parser.add_argument arg('--init', action='store_true', help='Initialize Berny optimizer.') arg('-f', '--format', choices=['xyz', 'aims'], default='xyz', help='Format of geometry') arg('-s', '--socket', nargs=2, metavar=('host', 'port'), help='Listen on given address') arg('paramfile', nargs='?', help='Optional optimization parameters as JSON') args = parser.parse_args() if args.init: init(args) elif args.socket: server(args) else: driver() PK!rk0G:G:berny/coords.py# This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. from __future__ import division from collections import OrderedDict from itertools import combinations, product import numpy as np from numpy import dot, pi from numpy.linalg import norm from . import Math from .species_data import get_property angstrom = 1/0.52917721092 #: class InternalCoord(object): def __init__(self, C=None): if C is not None: self.weak = sum( not C[self.idx[i], self.idx[i+1]] for i in range(len(self.idx)-1) ) def __eq__(self, other): self.idx == other.idx def __hash__(self): return hash(self.idx) def __repr__(self): args = list(map(str, self.idx)) if self.weak is not None: args.append('weak=' + str(self.weak)) return '{}({})'.format(self.__class__.__name__, ', '.join(args)) class Bond(InternalCoord): def __init__(self, i, j, **kwargs): if i > j: i, j = j, i self.i = i self.j = j self.idx = i, j InternalCoord.__init__(self, **kwargs) def hessian(self, rho): return 0.45*rho[self.i, self.j] def weight(self, rho, coords): return rho[self.i, self.j] def center(self, ijk): return np.round(ijk[[self.i, self.j]].sum(0)) def eval(self, coords, grad=False): v = (coords[self.i]-coords[self.j])*angstrom r = norm(v) if not grad: return r return r, [v/r, -v/r] class Angle(InternalCoord): def __init__(self, i, j, k, **kwargs): if i > k: i, j, k = k, j, i self.i = i self.j = j self.k = k self.idx = i, j, k InternalCoord.__init__(self, **kwargs) def hessian(self, rho): return 0.15*(rho[self.i, self.j]*rho[self.j, self.k]) def weight(self, rho, coords): f = 0.12 return np.sqrt(rho[self.i, self.j]*rho[self.j, self.k]) *\ (f+(1-f)*np.sin(self.eval(coords))) def center(self, ijk): return np.round(2*ijk[self.j]) def eval(self, coords, grad=False): v1 = (coords[self.i]-coords[self.j])*angstrom v2 = (coords[self.k]-coords[self.j])*angstrom dot_product = np.dot(v1, v2)/(norm(v1)*norm(v2)) if dot_product < -1: dot_product = -1 elif dot_product > 1: dot_product = 1 phi = np.arccos(dot_product) if not grad: return phi if abs(phi) > pi-1e-6: grad = [ (pi-phi)/(2*norm(v1)**2)*v1, (1/norm(v1)-1/norm(v2))*(pi-phi)/(2*norm(v1))*v1, (pi-phi)/(2*norm(v2)**2)*v2 ] else: grad = [ 1/np.tan(phi)*v1/norm(v1)**2-v2/(norm(v1)*norm(v2)*np.sin(phi)), (v1+v2)/(norm(v1)*norm(v2)*np.sin(phi)) - 1/np.tan(phi)*(v1/norm(v1)**2+v2/norm(v2)**2), 1/np.tan(phi)*v2/norm(v2)**2-v1/(norm(v1)*norm(v2)*np.sin(phi)) ] return phi, grad class Dihedral(InternalCoord): def __init__(self, i, j, k, l, weak=None, angles=None, C=None, **kwargs): if j > k: i, j, k, l = l, k, j, i self.i = i self.j = j self.k = k self.l = l self.idx = (i, j, k, l) self.weak = weak self.angles = angles InternalCoord.__init__(self, **kwargs) def hessian(self, rho): return 0.005*rho[self.i, self.j]*rho[self.j, self.k]*rho[self.k, self.l] def weight(self, rho, coords): f = 0.12 th1 = Angle(self.i, self.j, self.k).eval(coords) th2 = Angle(self.j, self.k, self.l).eval(coords) return (rho[self.i, self.j]*rho[self.j, self.k]*rho[self.k, self.l])**(1/3) * \ (f+(1-f)*np.sin(th1))*(f+(1-f)*np.sin(th2)) def center(self, ijk): return np.round(ijk[[self.j, self.k]].sum(0)) def eval(self, coords, grad=False): v1 = (coords[self.i]-coords[self.j])*angstrom v2 = (coords[self.l]-coords[self.k])*angstrom w = (coords[self.k]-coords[self.j])*angstrom ew = w/norm(w) a1 = v1-dot(v1, ew)*ew a2 = v2-dot(v2, ew)*ew sgn = np.sign(np.linalg.det(np.array([v2, v1, w]))) sgn = sgn or 1 dot_product = dot(a1, a2)/(norm(a1)*norm(a2)) if dot_product < -1: dot_product = -1 elif dot_product > 1: dot_product = 1 phi = np.arccos(dot_product)*sgn if not grad: return phi if abs(phi) > pi-1e-6: g = Math.cross(w, a1) g = g/norm(g) A = dot(v1, ew)/norm(w) B = dot(v2, ew)/norm(w) grad = [ g/(norm(g)*norm(a1)), -((1-A)/norm(a1)-B/norm(a2))*g, -((1+B)/norm(a2)+A/norm(a1))*g, g/(norm(g)*norm(a2)) ] elif abs(phi) < 1e-6: g = Math.cross(w, a1) g = g/norm(g) A = dot(v1, ew)/norm(w) B = dot(v2, ew)/norm(w) grad = [ g/(norm(g)*norm(a1)), -((1-A)/norm(a1)+B/norm(a2))*g, ((1+B)/norm(a2)-A/norm(a1))*g, -g/(norm(g)*norm(a2)) ] else: A = dot(v1, ew)/norm(w) B = dot(v2, ew)/norm(w) grad = [ 1/np.tan(phi)*a1/norm(a1)**2-a2/(norm(a1)*norm(a2)*np.sin(phi)), ((1-A)*a2-B*a1)/(norm(a1)*norm(a2)*np.sin(phi)) - 1/np.tan(phi)*((1-A)*a1/norm(a1)**2-B*a2/norm(a2)**2), ((1+B)*a1+A*a2)/(norm(a1)*norm(a2)*np.sin(phi)) - 1/np.tan(phi)*((1+B)*a2/norm(a2)**2+A*a1/norm(a1)**2), 1/np.tan(phi)*a2/norm(a2)**2-a1/(norm(a1)*norm(a2)*np.sin(phi)) ] return phi, grad def get_clusters(C): nonassigned = list(range(len(C))) clusters = [] while nonassigned: queue = set([nonassigned[0]]) clusters.append([]) while queue: node = queue.pop() clusters[-1].append(node) nonassigned.remove(node) queue.update(n for n in np.flatnonzero(C[node]) if n in nonassigned) C = np.zeros_like(C) for cluster in clusters: for i in cluster: C[i, cluster] = True return clusters, C class InternalCoords(object): def __init__(self, geom, allowed=None, dihedral=True, superweakdih=False): self._coords = [] n = len(geom) geom = geom.supercell() dist = geom.dist(geom) radii = np.array([get_property(sp, 'covalent_radius') for sp in geom.species]) bondmatrix = dist < 1.3*(radii[None, :]+radii[:, None]) self.fragments, C = get_clusters(bondmatrix) radii = np.array([get_property(sp, 'vdw_radius') for sp in geom.species]) shift = 0. C_total = C.copy() while not C_total.all(): bondmatrix |= ~C_total & (dist < radii[None, :]+radii[:, None]+shift) C_total = get_clusters(bondmatrix)[1] shift += 1. for i, j in combinations(range(len(geom)), 2): if bondmatrix[i, j]: bond = Bond(i, j, C=C) self.append(bond) for j in range(len(geom)): for i, k in combinations(np.flatnonzero(bondmatrix[j, :]), 2): ang = Angle(i, j, k, C=C) if ang.eval(geom.coords) > pi/4: self.append(ang) if dihedral: for bond in self.bonds: self.extend(get_dihedrals( [bond.i, bond.j], geom.coords, bondmatrix, C, superweak=superweakdih, )) if geom.lattice is not None: self._reduce(n) def append(self, coord): self._coords.append(coord) def extend(self, coords): self._coords.extend(coords) def __iter__(self): return self._coords.__iter__() def __len__(self): return len(self._coords) @property def bonds(self): return [c for c in self if isinstance(c, Bond)] @property def angles(self): return [c for c in self if isinstance(c, Angle)] @property def dihedrals(self): return [c for c in self if isinstance(c, Dihedral)] @property def dict(self): return OrderedDict([ ('bonds', self.bonds), ('angles', self.angles), ('dihedrals', self.dihedrals) ]) def __repr__(self): return "".format(', '.join( '{}: {}'.format(name, len(coords)) for name, coords in self.dict.items() )) def __str__(self): ncoords = sum(len(coords) for coords in self.dict.values()) s = 'Internal coordinates:\n' s += '* Number of fragments: {}\n'.format(len(self.fragments)) s += '* Number of internal coordinates: {}\n'.format(ncoords) for name, coords in self.dict.items(): for degree, adjective in [(0, 'strong'), (1, 'weak'), (2, 'superweak')]: n = len([None for c in coords if min(2, c.weak) == degree]) if n > 0: s += '* Number of {} {}: {}\n'.format(adjective, name, n) return s.rstrip() def eval_geom(self, geom, template=None): geom = geom.supercell() q = np.array([coord.eval(geom.coords) for coord in self]) if template is None: return q swapped = [] # dihedrals swapped by pi candidates = set() # potentially swapped angles for i, dih in enumerate(self): if not isinstance(dih, Dihedral): continue diff = q[i]-template[i] if abs(abs(diff)-2*pi) < pi/2: q[i] -= 2*pi*np.sign(diff) elif abs(abs(diff)-pi) < pi/2: q[i] -= pi*np.sign(diff) swapped.append(dih) candidates.update(dih.angles) for i, ang in enumerate(self): if not isinstance(ang, Angle) or ang not in candidates: continue # candidate angle was swapped if each dihedral that contains it was # either swapped or all its angles are candidates if all(dih in swapped or all(a in candidates for a in dih.angles) for dih in self.dihedrals if ang in dih.angles): q[i] = 2*pi-q[i] return q def _reduce(self, n): idxs = np.int64(np.floor(np.array(range(3**3*n))/n)) idxs, i = np.divmod(idxs, 3) idxs, j = np.divmod(idxs, 3) k = idxs % 3 ijk = np.vstack((i, j, k)).T-1 self._coords = [ coord for coord in self._coords if np.all(np.isin(coord.center(ijk), [0, -1])) ] idxs = set(i for coord in self._coords for i in coord.idx) self.fragments = [frag for frag in self.fragments if set(frag) & idxs] def hessian_guess(self, geom): geom = geom.supercell() rho = geom.rho() return np.diag([coord.hessian(rho) for coord in self]) def weights(self, geom): geom = geom.supercell() rho = geom.rho() return np.array([coord.weight(rho, geom.coords) for coord in self]) def B_matrix(self, geom): geom = geom.supercell() B = np.zeros((len(self), len(geom), 3)) for i, coord in enumerate(self): _, grads = coord.eval(geom.coords, grad=True) idx = [k % len(geom) for k in coord.idx] for j, grad in zip(idx, grads): B[i, j] += grad return B.reshape(len(self), 3*len(geom)) def update_geom(self, geom, q, dq, B_inv, log=lambda _: None): geom = geom.copy() thre = 1e-6 # target = CartIter(q=q+dq) # prev = CartIter(geom.coords, q, dq) for i in range(20): coords_new = geom.coords+B_inv.dot(dq).reshape(-1, 3)/angstrom dcart_rms = Math.rms(coords_new-geom.coords) geom.coords = coords_new q_new = self.eval_geom(geom, template=q) dq_rms = Math.rms(q_new-q) q, dq = q_new, dq-(q_new-q) if dcart_rms < thre: msg = 'Perfect transformation to cartesians in {} iterations' break if i == 0: keep_first = geom.copy(), q, dcart_rms, dq_rms else: msg = 'Transformation did not converge in {} iterations' geom, q, dcart_rms, dq_rms = keep_first log(msg.format(i+1)) log('* RMS(dcart): {:.3}, RMS(dq): {:.3}'.format(dcart_rms, dq_rms)) return q, geom def get_dihedrals(center, coords, bondmatrix, C, superweak=False): neigh_l = [n for n in np.flatnonzero(bondmatrix[center[0], :]) if n not in center] neigh_r = [n for n in np.flatnonzero(bondmatrix[center[-1], :]) if n not in center] angles_l = [Angle(i, center[0], center[1]).eval(coords) for i in neigh_l] angles_r = [Angle(center[-2], center[-1], j).eval(coords) for j in neigh_r] nonlinear_l = [n for n, ang in zip(neigh_l, angles_l) if ang < pi-1e-3 and ang >= 1e-3] nonlinear_r = [n for n, ang in zip(neigh_r, angles_r) if ang < pi-1e-3 and ang >= 1e-3] linear_l = [n for n, ang in zip(neigh_l, angles_l) if ang >= pi-1e-3 or ang < 1e-3] linear_r = [n for n, ang in zip(neigh_r, angles_r) if ang >= pi-1e-3 or ang < 1e-3] assert len(linear_l) <= 1 assert len(linear_r) <= 1 if center[0] < center[-1]: nweak = len(list( None for i in range(len(center)-1) if not C[center[i], center[i+1]] )) dihedrals = [] for nl, nr in product(nonlinear_l, nonlinear_r): if nl == nr: continue weak = nweak + \ (0 if C[nl, center[0]] else 1) + \ (0 if C[center[0], nr] else 1) if not superweak and weak > 1: continue dihedrals.append(Dihedral( nl, center[0], center[-1], nr, weak=weak, angles=( Angle(nl, center[0], center[1], C=C), Angle(nl, center[-2], center[-1], C=C) ) )) else: dihedrals = [] if len(center) > 3: pass elif linear_l and not linear_r: dihedrals.extend(get_dihedrals(linear_l + center, coords, bondmatrix, C)) elif linear_r and not linear_l: dihedrals.extend(get_dihedrals(center + linear_r, coords, bondmatrix, C)) return dihedrals PK!0T))berny/geomlib.py# This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. import os from io import StringIO from itertools import chain, groupby, product, repeat import numpy as np from numpy import pi from numpy.linalg import inv, norm from .species_data import get_property __version__ = '0.1.0' class Geometry(object): """ Represents a single molecule or a crystal. :param list species: list of element symbols :param list coords: list of atomic coordinates in angstroms (as 3-tuples) :param list lattice: list of lattice vectors (None for a moleucle) Iterating over a geometry yields 2-tuples of symbols and coordinates. ``len(geom)`` returns the number of atoms in a geometry. The class supports :py:func:`format` with the same available formats as :py:meth:`dump`. """ def __init__(self, species, coords, lattice=None): self.species = species self.coords = np.array(coords) self.lattice = np.array(lattice) if lattice is not None else None @classmethod def from_atoms(cls, atoms, lattice=None, unit=1.): """Alternative contructor. :param list atoms: list of 2-tuples with an elemnt symbol and a coordinate :param float unit: value to multiple atomic coordiantes with :param list lattice: list of lattice vectors (None for a moleucle) """ species = [sp for sp, _ in atoms] coords = [np.array(coord, dtype=float)*unit for _, coord in atoms] return cls(species, coords, lattice) def __repr__(self): s = repr(self.formula) if self.lattice is not None: s += ' in a lattice' return '<{} {}>'.format(self.__class__.__name__, s) def __iter__(self): for specie, coord in zip(self.species, self.coords): yield specie, coord def __len__(self): return len(self.species) @property def formula(self): """Chemical formula of the molecule or a unit cell.""" composition = sorted( (sp, len(list(g))) for sp, g in groupby(sorted(self.species)) ) return ''.join( '{}{}'.format(sp, n if n > 1 else '') for sp, n in composition ) def __format__(self, fmt): """Return the geometry represented as a string, delegates to :py:meth:`dump`.""" fp = StringIO() self.dump(fp, fmt) return fp.getvalue() dumps = __format__ def dump(self, f, fmt): """Saves the geometry into a file. :param file f: file object :param str fmt: geometry format, one of '', 'xyz', 'aims', 'mopac'. """ if fmt == '': f.write(repr(self)) elif fmt == 'xyz': f.write('{}\n'.format(len(self))) f.write('Formula: {}\n'.format(self.formula)) for specie, coord in self: f.write('{:>2} {}\n'.format( specie, ' '.join('{:15.8}'.format(x) for x in coord) )) elif fmt == 'aims': f.write('# Formula: {}\n'.format(self.formula)) for specie, coord in self: f.write('atom {} {:>2}\n'.format( ' '.join('{:15.8}'.format(x) for x in coord), specie )) elif fmt == 'mopac': f.write('* Formula: {}\n'.format(self.formula)) for specie, coord in self: f.write('{:>2} {}\n'.format( specie, ' '.join('{:15.8} 1'.format(x) for x in coord) )) else: raise ValueError("Unknown format: '{}'".format(fmt)) def copy(self): """Returns a copy of the geometry.""" return Geometry( list(self.species), self.coords.copy(), self.lattice.copy() if self.lattice is not None else None ) def write(self, filename): """ Writes the geometry into a file, delegates to :py:meth:`dump`. :param str filename: path that will be overwritten """ ext = os.path.splitext(filename)[1] if ext == '.xyz': fmt = 'xyz' elif ext == '.aims' or os.path.basename(filename) == 'geometry.in': fmt = 'aims' elif ext == '.mopac': fmt = 'mopac' else: raise ValueError('Unknown file extension') with open(filename, 'w') as f: self.dump(f, fmt) def super_circum(self, radius): """ Supercell dimensions such that the supercell circumsribes a sphere. :param float radius: circumscribed radius in angstroms Returns None when geometry is not a crystal. """ if self.lattice is None: return rec_lattice = 2*pi*inv(self.lattice.T) layer_sep = np.array( [sum(vec*rvec/norm(rvec)) for vec, rvec in zip(self.lattice, rec_lattice)] ) return np.array(np.ceil(radius/layer_sep+0.5), dtype=int) def supercell(self, ranges=((-1, 1), (-1, 1), (-1, 1)), cutoff=None): """ Creates a crystal supercell. :param list ranges: list of 2-tuples specifying the range of multiples of the unit-cell vectors :param float cutoff: if given, the ranges are determined such that the supercell contains a sphere with the radius qual to the cutoff Returns a copy of itself when geometry is not a crystal. """ if self.lattice is None: return self.copy() if cutoff: ranges = [(-r, r) for r in self.super_circum(cutoff)] latt_vectors = np.array([(0, 0, 0)] + [ sum(k*vec for k, vec in zip(shift, self.lattice)) for shift in product(*[range(a, b+1) for a, b in ranges]) if shift != (0, 0, 0) ]) species = list(chain.from_iterable(repeat(self.species, len(latt_vectors)))) coords = (self.coords[None, :, :]+latt_vectors[:, None, :]).reshape((-1, 3)) lattice = self.lattice*np.array([b-a for a, b in ranges])[:, None] return Geometry(species, coords, lattice) def dist_diff(self, other=None): r""" Calculate distances and vectors between atoms. :param Geometry other: calculate distances between two geometries if given or within a geometry if not Returns :math:`R_{ij}:=|\mathbf R_i-\mathbf R_j|` and :math:`R_{ij\alpha}:=(\mathbf R_i)_\alpha-(\mathbf R_j)_\alpha`. """ if other is None: other = self diff = self.coords[:, None, :]-other.coords[None, :, :] dist = np.sqrt(np.sum(diff**2, 2)) dist[np.diag_indices(len(self))] = np.inf return dist, diff def dist(self, other=None): """Returns the first element of :py:meth:`dist_diff`.""" return self.dist_diff(other)[0] def bondmatrix(self, scale=1.3): r""" Calculates the covalent connectedness matrix. :param float scale: threshold for accepting a distance as a covalent bond Returns :math:`b_{ij}:=R_{ij}<\text{scale}\times (R_i^\text{cov}+R_j^\text{cov})`. """ dist = self.dist(self) radii = np.array([get_property(sp, 'covalent_radius') for sp in self.species]) return dist < 1.3*(radii[None, :]+radii[:, None]) def rho(self): r""" Calculates a measure of covalentness. Returns :math:`\rho_{ij}:=\exp\big(-R_{ij}/(R_i^\text{cov}+R_j^\text{cov})\big)`. """ geom = self.supercell() dist = geom.dist(geom) radii = np.array([get_property(sp, 'covalent_radius') for sp in geom.species]) return np.exp(-dist/(radii[None, :]+radii[:, None])+1) @property def masses(self): """Returns an array of atomic masses.""" return np.array([get_property(sp, 'mass') for sp in self.species]) @property def cms(self): r"""Calculates the center of mass, :math:`\mathbf R_\text{CMS}`.""" masses = self.masses return np.sum(masses[:, None]*self.coords, 0)/masses.sum() @property def inertia(self): r"""Calculates the moment of inertia, :math:`I_{\alpha\beta}:= \sum_im_i\big(r_i^2\delta_{\alpha\beta}-(\mathbf r_i)_\alpha(\mathbf r_i)_\beta\big)` where :math:`\mathbf r_i=\mathbf R_i-\mathbf R_\text{CMS}`.""" coords_w = np.sqrt(self.masses)[:, None]*(self.coords-self.cms) A = np.array([np.diag(np.full(3, r)) for r in np.sum(coords_w**2, 1)]) B = coords_w[:, :, None]*coords_w[:, None, :] return np.sum(A-B, 0) def load(fp, fmt): """ Read a geometry from a file object. :param file fp: file object :param str fmt: the format of the geometry file, can be one of 'xyz', 'aims' Returns :py:class:`berny.Geometry`. """ if fmt == 'xyz': n = int(fp.readline()) fp.readline() species = [] coords = [] for _ in range(n): l = fp.readline().split() species.append(l[0]) coords.append([float(x) for x in l[1:4]]) return Geometry(species, coords) if fmt == 'aims': species = [] coords = [] lattice = [] while True: l = fp.readline() if l == '': break l = l.strip() if not l or l.startswith('#'): continue l = l.split() what = l[0] if what == 'atom': species.append(l[4]) coords.append([float(x) for x in l[1:4]]) elif what == 'lattice_vector': lattice.append([float(x) for x in l[1:4]]) if lattice: assert len(lattice) == 3 return Geometry(species, coords, lattice) else: return Geometry(species, coords) def loads(s, fmt): """ Read a geometry from a string, delegates to :py:func:`load`. :param str s: string with geometry """ fp = StringIO(s) return load(fp, fmt) def readfile(path, fmt=None): """ Read a geometry from a file path, delegates to :py:func:`load`. :param str path: path to a geometry file :param str fmt: if not given, the format is given from the file extension """ if not fmt: ext = os.path.splitext(path)[1] if ext == '.xyz': fmt = 'xyz' if ext == '.aims' or os.path.basename(path) == 'geometry.in': fmt = 'aims' with open(path) as f: return load(f, fmt) PK!\Vberny/optimizers.py# This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. __version__ = '0.2.0' def optimize(optimizer, solver): """ Optimize a geometry with respect to a solver. :param generator optimizer: Optimizer object with the same generator interface as :py:func:`berny.Berny` :param generator solver: unprimed generator that receives geometry as a 2-tuple of a list of 2-tuples of the atom symbol and coordinate (as a 3-tuple), and of a list of lattice vectors (or None if molecule), and yields the energy and gradients (as a N-by-3 matrix or (N+3)-by-3 matrix in case of a crystal geometry) Returns the optimized geometry. The function is equivalent to:: next(solver) for geom in optimizer: energy, gradients = solver.send((list(geom), geom.lattice)) optimizer.send((energy, gradients)) """ next(solver) for geom in optimizer: energy, gradients = solver.send((list(geom), geom.lattice)) optimizer.send((energy, gradients)) return geom PK!/w berny/solvers.py# Any copyright is dedicated to the Public Domain. # http://creativecommons.org/publicdomain/zero/1.0/ from __future__ import division import os import tempfile import subprocess import shutil import numpy as np from .coords import angstrom def MopacSolver(cmd='mopac', method='PM7', workdir=None): """ Wraps `MOPAC `_, which needs to be installed on the system. :param str cmd: MOPAC executable :param str method: model to calculate energy """ kcal, ev = 1/627.503, 1/27.2107 tmpdir = workdir or tempfile.mkdtemp() try: atoms, lattice = yield while True: mopac_input = '{} 1SCF GRADIENTS\n\n\n'.format(method) + '\n'.join( '{} {} 1 {} 1 {} 1'.format(el, *coord) for el, coord in atoms ) if lattice is not None: mopac_input += '\n' + '\n'.join( 'Tv {} 1 {} 1 {} 1'.format(*vec) for vec in lattice ) input_file = os.path.join(tmpdir, 'job.mop') with open(input_file, 'w') as f: f.write(mopac_input) subprocess.check_call([cmd, input_file]) with open(os.path.join(tmpdir, 'job.out')) as f: energy = float(next( l for l in f if 'TOTAL ENERGY' in l ).split()[3])*ev next(l for l in f if 'FINAL POINT AND DERIVATIVES' in l) next(f) next(f) gradients = np.array([ [float(next(f).split()[6])*kcal/angstrom for _ in range(3)] for _ in range(len(atoms)+(0 if lattice is None else 3)) ]) atoms, lattice = yield energy, gradients finally: if tmpdir != workdir: shutil.rmtree(tmpdir) def GenericSolver(f, *args, **kwargs): delta = kwargs.pop('delta', 1e-3) atoms, lattice = yield while True: energy = f(atoms, lattice, *args, **kwargs) coords = np.array([coord for _, coord in atoms]) gradients = np.zeros(coords.shape) for i_atom in range(coords.shape[0]): for i_xyz in range(3): ene = {} for step in [-2, -1, 1, 2]: coords_diff = coords.copy() coords_diff[i_atom, i_xyz] += step*delta atoms_diff = list(zip([sp for sp, _, in atoms], coords_diff)) ene[step] = f(atoms_diff, lattice, *args, **kwargs) gradients[i_atom, i_xyz] = _diff5(ene, delta) if lattice is not None: lattice_grads = np.zeros((3, 3)) for i_vec in range(3): for i_xyz in range(3): ene = {} for step in [-2, -1, 1, 2]: lattice_diff = lattice.copy() lattice_diff[i_vec, i_xyz] += step*delta ene[step] = f(atoms, lattice_diff, *args, **kwargs) lattice_grads[i_vec, i_xyz] = _diff5(ene, delta) gradients = np.vstack((gradients, lattice_grads)) atoms, lattice = yield energy, gradients/angstrom def _diff5(x, delta): return (1/12*x[-2]-2/3*x[-1]+2/3*x[1]-1/12*x[2])/delta PK!ʏberny/species-data.csv"number","name","symbol","covalent_radius","mass","vdw_radius" 1,"hydrogen","H",0.38,1.0079,1.6404493538 2,"helium","He",0.32,4.0026,1.4023196089 3,"lithium","Li",1.34,6.941,2.2013771973 4,"beryllium","Be",0.9,9.0122,2.2066689695 5,"boron","B",0.82,10.811,2.0584993504 6,"carbon","C",0.77,12.0107,1.8997461871 7,"nitrogen","N",0.75,14.0067,1.7674518844 8,"oxygen","O",0.73,15.9994,1.6880753028 9,"fluorine","F",0.71,18.9984,1.6086987211 10,"neon","Ne",0.69,20.1797,1.5399056837 11,"sodium","Na",1.54,22.9897,1.9738309967 12,"magnesium","Mg",1.3,24.305,2.2595866905 13,"aluminium","Al",1.18,26.9815,2.2913373232 14,"silicon","Si",1.11,28.0855,2.2225442858 15,"phosphorus","P",1.06,30.9738,2.1220006157 16,"sulfur","S",1.02,32.065,2.0426240341 17,"chlorine","Cl",0.99,35.453,1.9632474524 18,"argon","Ar",0.97,39.948,1.8785790987 19,"potassium","K",1.96,39.0983,1.9632474524 20,"calcium","Ca",1.74,40.078,2.4606740307 21,"scandium","Sc",1.44,44.9559,2.428923398 22,"titanium","Ti",1.36,47.867,2.3865892212 23,"vanadium","V",1.25,50.9415,2.3495468164 24,"chromium","Cr",1.27,51.9961,2.1114170715 25,"manganese","Mn",1.39,54.938,2.1008335273 26,"iron","Fe",1.25,55.845,2.2384196021 27,"cobalt","Co",1.26,58.9332,2.2119607416 28,"nickel","Ni",1.21,58.6934,2.0214569456 29,"copper","Cu",1.38,63.546,1.989706313 30,"zinc","Zn",1.31,65.39,2.1272923878 31,"gallium","Ga",1.26,69.723,2.2172525137 32,"germanium","Ge",1.22,72.64,2.2225442858 33,"arsenic","As",1.19,74.9216,2.1749183368 34,"selenium","Se",1.16,78.96,2.137875932 35,"bromine","Br",1.14,79.904,2.0796664388 36,"krypton","Kr",1.1,83.8,2.0214569456 37,"rubidium","Rb",2.11,85.4678,1.9685392245 38,"strontium","Sr",1.92,87.62,2.4024645375 39,"yttrium","Y",1.62,88.9059,2.5480411882 40,"zirconium","Zr",1.48,91.224,2.3971727654 41,"niobium","Nb",1.37,92.9064,2.241859254 42,"molybdenum","Mo",1.45,95.94,2.1690973875 43,"technetium","Tc",1.56,98,2.1569263116 44,"ruthenium","Ru",1.26,101.07,2.1142217107 45,"rhodium","Rh",1.35,102.9055,2.0902499831 46,"palladium","Pd",1.31,106.42,1.9367885919 47,"silver","Ag",1.53,107.8682,2.0214569456 48,"cadmium","Cd",1.48,112.411,2.1114170715 49,"indium","In",1.44,114.818,2.239467373 50,"tin","Sn",1.41,118.71,2.2770495385 51,"antimony","Sb",1.38,121.76,2.2627617538 52,"tellurium","Te",1.35,127.6,2.23312783 53,"iodine","I",1.33,126.9045,2.2066689695 54,"xenon","Xe",1.3,131.293,2.1590430205 55,"caesium","Cs",2.25,132.9055,2.0002898572 56,"barium","Ba",1.98,137.327,2.524175296 57,"lanthanum","La",1.69,138.9055,0 58,"cerium","Ce",,140.116,0 59,"praseodymium","Pr",,140.9077,0 60,"neodymium","Nd",,144.24,0 61,"promethium","Pm",,145,0 62,"samarium","Sm",,150.36,0 63,"europium","Eu",,151.964,0 64,"gadolinium","Gd",,157.25,0 65,"terbium","Tb",,158.9253,0 66,"dysprosium","Dy",,162.5,0 67,"holmium","Ho",,164.9303,0 68,"erbium","Er",,167.259,0 69,"thulium","Tm",,168.9342,0 70,"ytterbium","Yb",,173.04,0 71,"lutetium","Lu",1.6,174.967,0 72,"hafnium","Hf",1.5,178.49,2.2278360579 73,"tantalum","Ta",1.38,180.9479,2.1960854252 74,"tungsten","W",1.46,183.84,2.1590430205 75,"rhenium","Re",1.59,186.207,2.1272923878 76,"osmium","Os",1.28,190.23,2.0320404899 77,"iridium","Ir",1.37,192.217,2.1167088436 78,"platinum","Pt",1.28,195.078,2.0743746667 79,"gold","Au",1.44,196.9665,2.0426240341 80,"mercury","Hg",1.49,200.59,2.1061252994 81,"thallium","Tl",1.48,204.3833,2.0690828946 82,"lead","Pb",1.47,207.2,2.280753779 83,"bismuth","Bi",1.46,208.9804,2.2860455511 84,"polonium","Po",,209,2.1680390331 85,"astatine","At",,210,2.1537512484 86,"radon","Rn",1.45,222,2.2384196021 87,"francium","Fr",,223,0 88,"radium","Ra",,226,0 89,"actinium","Ac",,227,0 90,"thorium","Th",,232.0381,0 91,"protactinium","Pa",,231.0359,0 92,"uranium","U",,238.0289,0 PK!2Ԓ55berny/species_data.py# Any copyright is dedicated to the Public Domain. # http://creativecommons.org/publicdomain/zero/1.0/ import sys import csv from pkg_resources import resource_string def get_property(idx, name): if isinstance(idx, str): return species_data[idx][name] try: return next(row[name] for row in species_data if row['number'] == idx) except StopIteration: raise KeyError('No species with number "{}"'.format(idx)) def _get_species_data(): csv_lines = resource_string(__name__, 'species-data.csv').split(b'\n') if sys.version_info[0] > 2: csv_lines = [l.decode() for l in csv_lines] reader = csv.DictReader(csv_lines, quoting=csv.QUOTE_NONNUMERIC) species_data = {row['symbol']: row for row in reader} return species_data species_data = _get_species_data() PK!H#&((pyberny-0.4.1.dist-info/entry_points.txtN+I/N.,()JJ-ʫz9Vy\\PK!xFVAVApyberny-0.4.1.dist-info/LICENSEMozilla Public License Version 2.0 ================================== 1. Definitions -------------- 1.1. "Contributor" means each individual or legal entity that creates, contributes to the creation of, or owns Covered Software. 1.2. "Contributor Version" means the combination of the Contributions of others (if any) used by a Contributor and that particular Contributor's Contribution. 1.3. "Contribution" means Covered Software of a particular Contributor. 1.4. "Covered Software" means Source Code Form to which the initial Contributor has attached the notice in Exhibit A, the Executable Form of such Source Code Form, and Modifications of such Source Code Form, in each case including portions thereof. 1.5. "Incompatible With Secondary Licenses" means (a) that the initial Contributor has attached the notice described in Exhibit B to the Covered Software; or (b) that the Covered Software was made available under the terms of version 1.1 or earlier of the License, but not also under the terms of a Secondary License. 1.6. "Executable Form" means any form of the work other than Source Code Form. 1.7. "Larger Work" means a work that combines Covered Software with other material, in a separate file or files, that is not Covered Software. 1.8. "License" means this document. 1.9. "Licensable" means having the right to grant, to the maximum extent possible, whether at the time of the initial grant or subsequently, any and all of the rights conveyed by this License. 1.10. "Modifications" means any of the following: (a) any file in Source Code Form that results from an addition to, deletion from, or modification of the contents of Covered Software; or (b) any new file in Source Code Form that contains any Covered Software. 1.11. "Patent Claims" of a Contributor means any patent claim(s), including without limitation, method, process, and apparatus claims, in any patent Licensable by such Contributor that would be infringed, but for the grant of the License, by the making, using, selling, offering for sale, having made, import, or transfer of either its Contributions or its Contributor Version. 1.12. "Secondary License" means either the GNU General Public License, Version 2.0, the GNU Lesser General Public License, Version 2.1, the GNU Affero General Public License, Version 3.0, or any later versions of those licenses. 1.13. "Source Code Form" means the form of the work preferred for making modifications. 1.14. "You" (or "Your") means an individual or a legal entity exercising rights under this License. For legal entities, "You" includes any entity that controls, is controlled by, or is under common control with You. For purposes of this definition, "control" means (a) the power, direct or indirect, to cause the direction or management of such entity, whether by contract or otherwise, or (b) ownership of more than fifty percent (50%) of the outstanding shares or beneficial ownership of such entity. 2. License Grants and Conditions -------------------------------- 2.1. Grants Each Contributor hereby grants You a world-wide, royalty-free, non-exclusive license: (a) under intellectual property rights (other than patent or trademark) Licensable by such Contributor to use, reproduce, make available, modify, display, perform, distribute, and otherwise exploit its Contributions, either on an unmodified basis, with Modifications, or as part of a Larger Work; and (b) under Patent Claims of such Contributor to make, use, sell, offer for sale, have made, import, and otherwise transfer either its Contributions or its Contributor Version. 2.2. Effective Date The licenses granted in Section 2.1 with respect to any Contribution become effective for each Contribution on the date the Contributor first distributes such Contribution. 2.3. Limitations on Grant Scope The licenses granted in this Section 2 are the only rights granted under this License. No additional rights or licenses will be implied from the distribution or licensing of Covered Software under this License. Notwithstanding Section 2.1(b) above, no patent license is granted by a Contributor: (a) for any code that a Contributor has removed from Covered Software; or (b) for infringements caused by: (i) Your and any other third party's modifications of Covered Software, or (ii) the combination of its Contributions with other software (except as part of its Contributor Version); or (c) under Patent Claims infringed by Covered Software in the absence of its Contributions. This License does not grant any rights in the trademarks, service marks, or logos of any Contributor (except as may be necessary to comply with the notice requirements in Section 3.4). 2.4. Subsequent Licenses No Contributor makes additional grants as a result of Your choice to distribute the Covered Software under a subsequent version of this License (see Section 10.2) or under the terms of a Secondary License (if permitted under the terms of Section 3.3). 2.5. Representation Each Contributor represents that the Contributor believes its Contributions are its original creation(s) or it has sufficient rights to grant the rights to its Contributions conveyed by this License. 2.6. Fair Use This License is not intended to limit any rights You have under applicable copyright doctrines of fair use, fair dealing, or other equivalents. 2.7. Conditions Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted in Section 2.1. 3. Responsibilities ------------------- 3.1. Distribution of Source Form All distribution of Covered Software in Source Code Form, including any Modifications that You create or to which You contribute, must be under the terms of this License. You must inform recipients that the Source Code Form of the Covered Software is governed by the terms of this License, and how they can obtain a copy of this License. You may not attempt to alter or restrict the recipients' rights in the Source Code Form. 3.2. Distribution of Executable Form If You distribute Covered Software in Executable Form then: (a) such Covered Software must also be made available in Source Code Form, as described in Section 3.1, and You must inform recipients of the Executable Form how they can obtain a copy of such Source Code Form by reasonable means in a timely manner, at a charge no more than the cost of distribution to the recipient; and (b) You may distribute such Executable Form under the terms of this License, or sublicense it under different terms, provided that the license for the Executable Form does not attempt to limit or alter the recipients' rights in the Source Code Form under this License. 3.3. Distribution of a Larger Work You may create and distribute a Larger Work under terms of Your choice, provided that You also comply with the requirements of this License for the Covered Software. If the Larger Work is a combination of Covered Software with a work governed by one or more Secondary Licenses, and the Covered Software is not Incompatible With Secondary Licenses, this License permits You to additionally distribute such Covered Software under the terms of such Secondary License(s), so that the recipient of the Larger Work may, at their option, further distribute the Covered Software under the terms of either this License or such Secondary License(s). 3.4. Notices You may not remove or alter the substance of any license notices (including copyright notices, patent notices, disclaimers of warranty, or limitations of liability) contained within the Source Code Form of the Covered Software, except that You may alter any license notices to the extent required to remedy known factual inaccuracies. 3.5. Application of Additional Terms You may choose to offer, and to charge a fee for, warranty, support, indemnity or liability obligations to one or more recipients of Covered Software. However, You may do so only on Your own behalf, and not on behalf of any Contributor. You must make it absolutely clear that any such warranty, support, indemnity, or liability obligation is offered by You alone, and You hereby agree to indemnify every Contributor for any liability incurred by such Contributor as a result of warranty, support, indemnity or liability terms You offer. You may include additional disclaimers of warranty and limitations of liability specific to any jurisdiction. 4. Inability to Comply Due to Statute or Regulation --------------------------------------------------- If it is impossible for You to comply with any of the terms of this License with respect to some or all of the Covered Software due to statute, judicial order, or regulation then You must: (a) comply with the terms of this License to the maximum extent possible; and (b) describe the limitations and the code they affect. Such description must be placed in a text file included with all distributions of the Covered Software under this License. Except to the extent prohibited by statute or regulation, such description must be sufficiently detailed for a recipient of ordinary skill to be able to understand it. 5. Termination -------------- 5.1. The rights granted under this License will terminate automatically if You fail to comply with any of its terms. However, if You become compliant, then the rights granted under this License from a particular Contributor are reinstated (a) provisionally, unless and until such Contributor explicitly and finally terminates Your grants, and (b) on an ongoing basis, if such Contributor fails to notify You of the non-compliance by some reasonable means prior to 60 days after You have come back into compliance. Moreover, Your grants from a particular Contributor are reinstated on an ongoing basis if such Contributor notifies You of the non-compliance by some reasonable means, this is the first time You have received notice of non-compliance with this License from such Contributor, and You become compliant prior to 30 days after Your receipt of the notice. 5.2. If You initiate litigation against any entity by asserting a patent infringement claim (excluding declaratory judgment actions, counter-claims, and cross-claims) alleging that a Contributor Version directly or indirectly infringes any patent, then the rights granted to You by any and all Contributors for the Covered Software under Section 2.1 of this License shall terminate. 5.3. In the event of termination under Sections 5.1 or 5.2 above, all end user license agreements (excluding distributors and resellers) which have been validly granted by You or Your distributors under this License prior to termination shall survive termination. ************************************************************************ * * * 6. Disclaimer of Warranty * * ------------------------- * * * * Covered Software is provided under this License on an "as is" * * basis, without warranty of any kind, either expressed, implied, or * * statutory, including, without limitation, warranties that the * * Covered Software is free of defects, merchantable, fit for a * * particular purpose or non-infringing. The entire risk as to the * * quality and performance of the Covered Software is with You. * * Should any Covered Software prove defective in any respect, You * * (not any Contributor) assume the cost of any necessary servicing, * * repair, or correction. This disclaimer of warranty constitutes an * * essential part of this License. No use of any Covered Software is * * authorized under this License except under this disclaimer. * * * ************************************************************************ ************************************************************************ * * * 7. Limitation of Liability * * -------------------------- * * * * Under no circumstances and under no legal theory, whether tort * * (including negligence), contract, or otherwise, shall any * * Contributor, or anyone who distributes Covered Software as * * permitted above, be liable to You for any direct, indirect, * * special, incidental, or consequential damages of any character * * including, without limitation, damages for lost profits, loss of * * goodwill, work stoppage, computer failure or malfunction, or any * * and all other commercial damages or losses, even if such party * * shall have been informed of the possibility of such damages. This * * limitation of liability shall not apply to liability for death or * * personal injury resulting from such party's negligence to the * * extent applicable law prohibits such limitation. Some * * jurisdictions do not allow the exclusion or limitation of * * incidental or consequential damages, so this exclusion and * * limitation may not apply to You. * * * ************************************************************************ 8. Litigation ------------- Any litigation relating to this License may be brought only in the courts of a jurisdiction where the defendant maintains its principal place of business and such litigation shall be governed by laws of that jurisdiction, without reference to its conflict-of-law provisions. Nothing in this Section shall prevent a party's ability to bring cross-claims or counter-claims. 9. Miscellaneous ---------------- This License represents the complete agreement concerning the subject matter hereof. If any provision of this License is held to be unenforceable, such provision shall be reformed only to the extent necessary to make it enforceable. Any law or regulation which provides that the language of a contract shall be construed against the drafter shall not be used to construe this License against a Contributor. 10. Versions of the License --------------------------- 10.1. New Versions Mozilla Foundation is the license steward. Except as provided in Section 10.3, no one other than the license steward has the right to modify or publish new versions of this License. Each version will be given a distinguishing version number. 10.2. Effect of New Versions You may distribute the Covered Software under the terms of the version of the License under which You originally received the Covered Software, or under the terms of any subsequent version published by the license steward. 10.3. Modified Versions If you create software not governed by this License, and you want to create a new license for such software, you may create and use a modified version of this License if you rename the license and remove any references to the name of the license steward (except to note that such modified license differs from this License). 10.4. Distributing Source Code Form that is Incompatible With Secondary Licenses If You choose to distribute Source Code Form that is Incompatible With Secondary Licenses under the terms of this version of the License, the notice described in Exhibit B of this License must be attached. Exhibit A - Source Code Form License Notice ------------------------------------------- This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of the MPL was not distributed with this file, You can obtain one at http://mozilla.org/MPL/2.0/. If it is not possible or desirable to put the notice in a particular file, then You may include the notice in a location (such as a LICENSE file in a relevant directory) where a recipient would be likely to look for such a notice. You may add additional accurate notices of copyright ownership. Exhibit B - "Incompatible With Secondary Licenses" Notice --------------------------------------------------------- This Source Code Form is "Incompatible With Secondary Licenses", as defined by the Mozilla Public License, v. 2.0. PK!HfWXpyberny-0.4.1.dist-info/WHEEL A н#f@Z|Jt~7m7 G9oF04hAf"Rhzi;0k:؂{ZO, JPK!HT" pyberny-0.4.1.dist-info/METADATAVmO#7_aGNYox+P);*CUP%c?c&P("<3c9<Zetϼv Vd4b[x5mFtMmu676^,kH /`qYҗ͘ ]tt*(ȳdAKm3Wؚ+%PsYe4/߸*E&|oD6ΐm ֧`} wNJ@c5VYFrDyP9%(x)0܊rp_.O1VO |1`P9C[]X^wUGd :6ų%Ҍ?-͘J)f'qVǪ +!J%&{,lWc?x3*\ŋ9?%0ng A;[$6m`u` ! j{:Z_uꡳә}^P@ w{h@iFhBebh55MX^9׶ݑ +MI`5r}Y7qCq#ϵ J UDyPenR 2i[AgJ Rч>ocA^Q.\XJJH7n7| ?&O^+, 0z﨓؎_I/@C/h1 }ig-4?)$G,ߔ6͚;>Б1 dP^u ; ?AFH)aX+ooFL: 7a0o~ ]&-99Qh,p*2]lRK,Bx翙3,X2zZATlmfezsB0k9[DkC r1h,x1ͱ1`+YwkML)p ٦ q T#EQL̀v d+ 2{=r̚04KEG'aͺP'yUQ%_0V2:Tښӝ_}v'xD$ðo( #O QMSYSMxѤMx>sؓomGCBb$#*PK!#  berny/Logger.pyPK!w( 9berny/Math.pyPK!&pLberny/__init__.pyPK!vYqXX !berny/_py2.pyPK!T2%%berny/berny.pyPK!JPB B ;berny/cli.pyPK!rk0G:G:WFberny/coords.pyPK!0T))ˀberny/geomlib.pyPK!\VӪberny/optimizers.pyPK!/w berny/solvers.pyPK!ʏberny/species-data.csvPK!2Ԓ55{berny/species_data.pyPK!H#&((pyberny-0.4.1.dist-info/entry_points.txtPK!xFVAVAOpyberny-0.4.1.dist-info/LICENSEPK!HfWXpyberny-0.4.1.dist-info/WHEELPK!HT" tpyberny-0.4.1.dist-info/METADATAPK!H#] qpyberny-0.4.1.dist-info/RECORDPKq