From 4bdbbc6a0413518999dfe19c2e5acda946c0d0af Mon Sep 17 00:00:00 2001 From: Amal Kanta Giri Date: Mon, 19 Sep 2022 16:23:58 +0200 Subject: [PATCH] Contact angle observable reimplemented * implementation of instantaneous contact angles and other geometrical parameters * implementation of ellipsoidal linear and nonlinear fitting, with resolution of contact line location and contact angle function * new linear least squares and nonlinear fitting of ellipse, with RMSD tested to be lower than the fit to a circumference * ellipticity constraint imposed for linear and nonlinear fit * new procedure to remove automatically the liquid/solid interface atoms and obtain the liquid/vapur interface * more consistent user interface --- pytim/observables/__init__.py | 2 +- pytim/observables/contactangle.py | 1496 ++++++++++++++++++++--------- 2 files changed, 1043 insertions(+), 455 deletions(-) diff --git a/pytim/observables/__init__.py b/pytim/observables/__init__.py index 8661f414..1785e0a9 100644 --- a/pytim/observables/__init__.py +++ b/pytim/observables/__init__.py @@ -7,7 +7,7 @@ from MDAnalysis.core.groups import Atom, AtomGroup, Residue, ResidueGroup from .observable import Observable -from .contactangle import ContactAngleGitim, ContactAngleGibbs +from .contactangle import ContactAngle from .basic_observables import Position, RelativePosition, Velocity, Force from .basic_observables import Number, Mass, Charge, NumberOfResidues from .basic_observables import Distance diff --git a/pytim/observables/contactangle.py b/pytim/observables/contactangle.py index 90a6c8e0..15c2560e 100644 --- a/pytim/observables/contactangle.py +++ b/pytim/observables/contactangle.py @@ -3,524 +3,1112 @@ """ Module: ContactAngle ==================== """ - +import pytim import MDAnalysis as mda import numpy as np from scipy.spatial import cKDTree +import scipy.linalg +import warnings +from numpy.linalg import solve +from scipy.optimize import curve_fit +from scipy.optimize import minimize class ContactAngle(object): - """ Base class for contact angle calculation - """ + """ ContactAngle class implementation that uses interfacial atoms to compute + droplet profiles and contact angles using different approaches. + + :param PYTIM interface : Compute the contact angle for this interface + :param AtomGroup droplet : Atom group representative of the droplet + :param AtomGroup substrate : Atom group representative of the substrate + :param int periodic : direction along which the system is periodic. Default: None, not periodic + If None, the code performs best fit to ellipsoids, otherwise, to ellipses + If not None, selects the direction of the axis of macroscopic translational + invariance: 0:x, 1:y, 2:z + :param float hcut : don't consider contributions from atoms closer than this to the substrate + (used to disregard the microscopic contact angle) + :param float hcut_upper : don't consider contributions from atoms above this distance from the substrate + default: None + :param int bins : bins used for sampling the profile + :param int or array_like removeCOM : remove the COM motion along this direction(s). Default: None, does not remove COM motion + + Example: + + >>> import MDAnalysis as mda + >>> import numpy as np + >>> import pytim + >>> from pytim import observables + >>> from pytim.datafiles import * + >>> + >>> u = mda.Universe(WATER_DROPLET_CYLINDRICAL_GRO,WATER_DROPLET_CYLINDRICAL_XTC) + >>> droplet = u.select_atoms("name OW") + >>> substrate = u.select_atoms("name C") + >>> inter = pytim.GITIM(universe=u,group=droplet, molecular=False,alpha=2.5,cluster_cut=3.4, biggest_cluster_only=True) + >>> + >>> # Contact angle calculation using interfacial atoms, angular bining for a cylindrical droplet + >>> # with periodicity along the y (1) axis + >>> + >>> CA = observables.ContactAngle(inter, substrate, periodic=1,bins=397,removeCOM=0,hcut=5) + >>> for ts in u.trajectory[::]: + ... CA.sample() + >>> # Instantaneous contact angle (last frame) by fitting a circle... + >>> np.round(CA.contact_angle,2) + 90.58 + + >>> + >>> # ... and using an elliptical fit: + >>> left, right = CA.contact_angles + >>> # left angle + >>> np.round(np.abs(left),2) + 79.95 + + >>> # right angle + >>> np.round(right,2) + 83.84 + + >>> # Contact angles from the averaged binned statistics of + >>> # surface atoms' radial distance as a function of the azimuthal angle + >>> list(np.round(CA.mean_contact_angles,2)) + [96.2, 100.68] + + """ + + class Histogram(object): + def __init__(self): + self.x, self.y, self.h = np.array([]), np.array([]), np.array([]) + + histo = Histogram() + + def __init__(self, inter, substrate, periodic=None, hcut=0.0, hcut_upper=None, + contact_cut=0.0, bins=100, removeCOM=None, store=False): + self._sampled = False + self.inter, self.substrate = inter, substrate + self.universe, self.trajectory = self.inter.universe, self.inter.universe.trajectory + self.hcut, self.hcut_upper = hcut, hcut_upper + self.nframes, self.periodic = 0, periodic + self.bins, self.removeCOM = bins, removeCOM + self.contact_cut, self.store = contact_cut, store + self.x, self.y, self.h = np.array([]), np.array([]), np.array( + []) # in 2d, x and y are radius and elevation + if self.contact_cut == 'auto': + tree = cKDTree(substrate.positions, + boxsize=substrate.universe.dimensions[:3]) + d, _ = tree.query(inter.atoms.positions, k=1) + hist, bins = np.histogram(d, int(np.max(d)/0.5)) + # check the location of the maximum with resolution 0.5 Angstrom + vmax, imax = np.max(hist), np.argmax(hist) + # find where the distribution drops below 5% of the maximum; this could happen before and after the maximum + idx = np.where(hist <= 0.05*vmax)[0] + # select the first bin that satisfies the above condition and is after the maximum + self.contact_cut = bins[idx[idx > imax][0]] + + @property + def contact_angle(self): + """ + The contact angle from the best-fititng circumference or sphere + computed using the current frame + """ + if self.periodic is None: + raise NotImplementedError( + 'fit_sphere not implemented, use contact_angles instead of contact_angle') + else: + c, cp, theta, _, _ = self.fit_circle() + self._polynomial_coefficients = c + self._canonical_form = cp + return theta + + @property + def mean_contact_angle(self): + """ + The contact angle from the best-fititng circumference + computed using the location of the interface averaged + over the sampled frames + """ + if self.periodic is None: + raise NotImplementedError( + 'fit_sphere not implemented, use contact_angles instead of contact_angle') + else: + c, cp, theta, _, _ = self.fit_circle(use='histogram') + self._polynomial_coefficients = c + self._canonical_form = cp + return theta - # todo pass surface normal and adapt algorithms. - def __init__(self, droplet, substrate, zcut=5.0, debug=False): + @property + def contact_angles(self): + """ + The contact angles from the best-fititng ellipse or + ellipsoid computed using the current frame """ - Computes radial profiles and contact angles using different approaches - - :param AtomGroup fluid : Atom group representative of the droplet - :param AtomGroup substrate : Atom group representative of the substrate - :param float zcut : don't consider contributions from atoms nearer than this to the substrate - - Example: - - >>> import MDAnalysis as mda - >>> import numpy as np - >>> import pytim - >>> from pytim import observables - >>> from pytim.datafiles import * - >>> - >>> u = mda.Universe(WATER_DROPLET_CYLINDRICAL_GRO,WATER_DROPLET_CYLINDRICAL_XTC) - >>> droplet = u.select_atoms("name OW") - >>> substrate = u.select_atoms("name C") - >>> - >>> CA = observables.ContactAngleGitim(droplet, substrate) - >>> - >>> inter = pytim.GITIM(universe=u,group=droplet, molecular=False,alpha=2.5,cluster_cut=3.4, biggest_cluster_only=True) - >>> - >>> for ts in u.trajectory[::]: - ... CA.sample(inter) - >>> - >>> np.round(CA.contact_angle,2) - 79.52 + if self.periodic is None: + c, cp, theta, _ = self.fit_ellipsoid() + else: + c, cp, theta, _ = self.fit_ellipse() + self._polynomial_coefficients = c + self._canonical_form = cp + return theta + @property + def mean_contact_angles(self): + """ + The contact angles from the best-fititng ellipse or + ellipsoid computed using the location of the interface + averaged over the sampled frames """ + if self.periodic is None: + c, cp, theta, _ = self.fit_ellipsoid(use='histogram') + else: + c, cp, theta, _ = self.fit_ellipse(use='histogram') + self._polynomial_coefficients = c + self._canonical_form = cp + return theta - self.droplet = droplet - self.substrate = substrate - self.universe = self.droplet.universe - self.trajectory = self.droplet.universe.trajectory - self.zcut = zcut - self.nframes = 0 - self.debug = debug - masses = self.universe.atoms.masses.copy() - masses[self.universe.atoms.masses == 0.0] = 40.0 - self.universe.atoms.masses = masses.copy() - - def remove_com(self, group): - if id(self.droplet.universe) != id(group.universe): - raise RuntimeError('Universes are not the same') - com = self.droplet.center_of_mass() - zoff = np.max(self.substrate.atoms.positions[:, 2]) - com[2] = zoff + @property + def polynomial_coefficients(self): + try: + return self._polynomial_coefficients + except: + raise RuntimeError('No fit has been performed yet') + + @property + def canonical_form(self): + try: + return self._canonical_form + except: + raise RuntimeError('No fit has been performed yet') + + def _shifted_pos(self, group): + com = group.center_of_mass() + hoff = np.max(self.substrate.atoms.positions[:, 2]) + com[2] = hoff + # removes the com along x,y and uses the topmost substrate atom to + # define the location of the origin along the vertical axis pos = group.positions - com - pos = pos[pos[:, 2] > self.zcut] + # further, we won't be using any coordinate that is below a certain + # distance from the substrate surface + pos = pos[pos[:, 2] > self.hcut] + if self.hcut_upper is not None: + pos = pos[pos[:, 2] < self.hcut_upper] return pos + def _remove_COM(self, group, inter, alpha, box): + def group_size(g, direction): + p = g.atoms.positions[:, direction] + return np.abs(np.max(p) - np.min(p)) + + if self.removeCOM is not None: + RC = self.removeCOM + if type(RC) == int: + RC = [RC] + for axis in RC: + while (group_size(group, axis) > box[axis] - 4 * alpha): + pos = group.universe.atoms.positions.copy() + pos[:, axis] += alpha + group.universe.atoms.positions = pos.copy() + pos = group.universe.atoms.pack_into_box() + group.universe.atoms.positions = pos.copy() + com = self.inter.atoms.center_of_mass()[axis] + pos = group.universe.atoms.positions.copy() + pos[:, axis] -= com + group.universe.atoms.positions = pos.copy() + + pos[:, axis] += box[axis]/2 + group.universe.atoms.positions = pos.copy() + + pos = group.universe.atoms.pack_into_box() + group.universe.atoms.positions = pos.copy() + return group.universe.atoms.positions + + def sample(self): + """ samples the profile of a droplet, stores + the current atomic coordinates of the liquid surface + not in contact with the substrate in the reference + frame of the droplet, and optionally the coordinates + along the whole trajectory. + """ + if self.contact_cut > 0: + tree_substrate = cKDTree( + self.substrate.atoms.positions, boxsize=self.substrate.universe.dimensions[:3]) + dists, _ = tree_substrate.query(self.inter.atoms.positions, k=1) + self.liquid_vapor = self.inter.atoms[dists > self.contact_cut] + else: + self.liquid_vapor = self.inter.atoms + + box = self.substrate.universe.dimensions[:3] + if self.store: + try: + _ = self._x + except: + self._x, self._y, self._h = [], [], [] + + self.nframes += 1 + + self._remove_COM(self.inter.atoms, self.liquid_vapor, + self.inter.alpha, box) + self.maxr = np.max(self.substrate.universe.dimensions[:2]) + self.maxh = self.substrate.universe.dimensions[2] + + self.x, self.y, self.h = self.cartesian_coordinates( + self.liquid_vapor.atoms) + if self.store: + self._x += list(self.x) + self._h += list(self.h) + if self.periodic is None: + self._y += list(self.y) + + if self.periodic is None: + r, theta, h, phi = self.spherical_coordinates( + self.liquid_vapor.atoms) + self.phi = phi.copy() + else: + r, theta, h, z = self.cylindrical_coordinates( + self.liquid_vapor.atoms) + + self.r = r.copy() + self.theta = theta.copy() + self.h = h.copy() + + if self.periodic is None: + pass # TODO implement the 2d histogram + else: + self.sample_theta_R(theta, r, h) + @staticmethod - def arc(r, center, rad): - return center+np.sqrt(rad**2-r**2) + def rmsd_ellipse(p, x, N, check_coeffs=True): + cp = ContactAngle._ellipse_general_to_canonical(p, check_coeffs) + xx, yy = ContactAngle.ellipse(cp, N) + pos = np.vstack([xx, yy]).T + cond = np.logical_and(yy >= np.min(x[:, 1]), yy <= np.max(x[:, 1])) + pos = pos[cond] + return np.sqrt(np.mean(cKDTree(pos).query(x)[0]**2)) @staticmethod - def fit_arc_(z, r, rmin=None, rmax=None, use_points=False, p0=None): - """ fit an arc through the profile z(r) sampled by the class + def rmsd_ellipse_penalty(p, x, N, check_coeffs=True): + rmsd = ContactAngle.rmsd_ellipse(p, x, N, check_coeffs) + penalty = (4*p[0]*p[2]-p[1]**2-1)**2 + return rmsd + penalty - :param float rmin : minimum radius used for the fit (default: 0) - :param float rmax : maximum radius used for the fit (default: maximum from dataset) - :param bool use_points : False: use the profile estimate, True: use the surface atoms positions (only) + @staticmethod + def rmsd_ellipsoid(p, x, N, check_coeffs=True): + """ RMSD between the points x and the ellipsoid defined by the general + parameters p of the associated polynomial. + + :param p list : general coefficients [a,b,c,f,g,h,p,q,r,d] + :param x ndarray : points coordinates as a (*,3)-ndarray + :param N int : number of points on the ellipsoid that are + generated and used to compute the rmsd + :param check_coeffs bool : if true, perform additional checks """ - def arc(r, center, R): - return center+np.sqrt(R**2-r**2) - - from scipy.optimize import curve_fit - - z = np.asarray(z) - r = np.asarray(r) - - con = np.logical_and(np.isfinite(z), np.isfinite(r)) - z, r = z[con], r[con] - - if rmax is None: - rmax = np.nanmax(r) - rmin = np.nanmin(r) - cond = np.logical_and(r > rmin, r < rmax) - if p0 is None: - p0 = (-20., 110.) - popt, pcov = curve_fit(arc, r[cond], z[cond], p0=p0) - center, rad = popt - #print("center=",center, "radius=",rad) - rad = np.abs(rad) - base_radius = np.sqrt(rad**2-center**2) - # old from marcello costheta = -np.sin(center/rad) - costheta = -center/rad - return (rad, base_radius, costheta, center) + cp = ContactAngle._ellipsoid_general_to_affine(p, check_coeffs) + xx, yy, zz = ContactAngle.ellipsoid(cp, N) + pos = np.vstack([xx, yy, zz]).T + cond = np.logical_and(zz >= np.min(x[:, 2]), zz <= np.max(x[:, 2])) + pos = pos[cond] + return np.sqrt(np.mean(cKDTree(pos).query(x)[0]**2)) @staticmethod - def fit_ellipse_(r, z, yy=None): - ZZ = np.asarray(z) - RR = np.asarray(r) - idx1 = np.isnan(RR) - idx2 = np.isnan(ZZ) - idx = np.logical_and(~idx1, ~idx2) # nan removed - R_, Z_ = RR[idx], ZZ[idx] - x1, y1 = R_, Z_ - X = x1[:, np.newaxis] - Y = y1[:, np.newaxis] - # Formulate and solve the least squares problem ||Ax - b ||^2 - A = np.hstack([X**2, X * Y, Y**2, X, Y]) - b = np.ones_like(X) - a = np.linalg.lstsq(A, b)[0].squeeze() - # Print the equation of the ellipse in standard form - print('Ellipse equ: {0:.3}x^2 + {1:.3}xy+{2:.3}y^2+{3:.3}x+{4:.3}y = 0'.format( - a[0], a[1], a[2], a[3], a[4])) - if yy is None: - yy = 0 # positon of substrate - bc = a[1]*yy+a[3] - cc = a[2]*yy**2+a[4]*yy-1 - x1 = (-bc+np.sqrt(bc**2-4*a[0]*cc))/(2*a[0]) - x2 = (-bc-np.sqrt(bc**2-4*a[0]*cc))/(2*a[0]) - m1 = -(2*a[0]*x1+a[3]+a[1]*yy)/(a[1]*x1+a[4]+2*a[2]*yy) - m2 = -(2*a[0]*x2+a[3]+a[1]*yy)/(a[1]*x2+a[4]+2*a[2]*yy) - theta1 = np.arctan(m1)*180/np.pi - theta2 = np.arctan(m2)*180/np.pi - return a, theta1, theta2 + def rmsd_ellipsoid_penalty(p, x, N, check_coeffs=True): + rmsd = ContactAngle.rmsd_ellipsoid(p, x, N, check_coeffs) + violation = (4*p[0]*p[2]-p[1]**2)**2 + penalty = 0 if 4*p[0]*p[2]-p[1]**2 > 0 else violation + return rmsd + penalty + + @staticmethod + def rmsd_circle(p, x): + R, x0, y0 = p + d = np.linalg.norm(np.array([x0, y0])-x, axis=1) - R + return np.sqrt(np.mean(d**2)) @staticmethod - def arc_ellipse(x, rmin=None, rmax=None, bins=None): - if rmax is None: - rmax, rmin, bins = 80, -80, 3200 - #x_coord = np.linspace(rmin,rmax,bins) - #y_coord = np.linspace(zmin,zmax,bins) - #X_coord, Y_coord = np.meshgrid(x_coord, y_coord) - #Z_coord = x[0] * X_coord ** 2 + x[1] * X_coord * Y_coord + x[2] * Y_coord**2 + x[3] * X_coord + x[4] * Y_coord - # rmin,rmax,bins=int(rmin),int(rmax),int(bins) - val = np.linspace(rmin, rmax, bins) - bb = x[1]*val+x[4] - # aa=x[2] - cc = x[0]*val*val+x[3]*val-1 - yyP = (-bb+np.sqrt(bb*bb - 4*x[2]*cc))/(2*x[2]) - yyN = (-bb-np.sqrt(bb*bb - 4*x[2]*cc))/(2*x[2]) - yy = list(yyP)+list(yyN) - xx = list(val)+list(val) - xx = np.asarray(xx) - yy = np.asarray(yy) - idy = ~np.isnan(yy) - # return X_coord,Y_coord,Z_coord,xx[idy],yy[idy] - return xx[idy], yy[idy] - - def fit_arc(self, rmin=0, rmax=None, use_points=None, p0=None): - """ fit an arc through the profile z(r) sampled by the class - - :param float rmin : minimum radius used for the fit (default: 0) - :param float rmax : maximum radius used for the fit (default: maximum from dataset) + def ellipse(parmsc, npts=100, tmin=0., tmax=2.*np.pi): + """ Return npts points on the ellipse described by the canonical parameters + x0, y0, ap, bp, e, phi for values of the paramter between tmin and tmax. + + :param dict parmsc : dictionary with keys: x0,y0,a,b,phi + :param float tmin : minimum value of the parameter + :param float tmax : maximum value of the parameter + :param int npts : number of points to use + return: + :tuple : (x,y) of coordinates as numpy arrays """ - try: - z, r = self.z, self.r + t = np.linspace(tmin, tmax, npts) + x = parmsc['x0'] + parmsc['a'] * np.cos(t) * np.cos( + parmsc['phi']) - parmsc['b'] * np.sin(t) * np.sin(parmsc['phi']) + y = parmsc['y0'] + parmsc['a'] * np.cos(t) * np.sin( + parmsc['phi']) + parmsc['b'] * np.sin(t) * np.cos(parmsc['phi']) + return x, y - except AttributeError: - raise RuntimeError( - 'something is wrong: no surface atoms or not gibbs surface present') + @staticmethod + def circle(parmsc, npts=100, tmin=0., tmax=2.*np.pi): + """ Return npts points on the circle described by the canonical parameters + R, x0, y0 for values of the paramter between tmin and tmax. - return self.fit_arc_(z, r, rmin=rmin, rmax=rmax, p0=p0) + :param dict parmsc : dictionary with keys: R, x0, y0 + :param float tmin : minimum value of the parameter + :param float tmax : maximum value of the parameter + :param int npts : number of points to use - def fit_arc_ellipse(self): - try: - z, r = self.z, self.r - except AttributeError: - raise RuntimeError( - 'something is wrong: no surface atoms or not gibbs surface present') + return: + :tuple : (x,y) of coordinates as numpy arrays + """ + t = np.linspace(tmin, tmax, npts) + x = parmsc['x0'] + parmsc['R'] * np.cos(t) + y = parmsc['y0'] + parmsc['R'] * np.sin(t) + return x, y + + @staticmethod + def ellipsoid(parmsc, npts=1000): + """ Return npts points on the ellipsoid described by the affine parameters + T, v - return self.fit_ellipse_(r, z, yy=None) + :param dict parmsc : dictionary with keys: T (3x3 matrix), v (3x1 vector) + :param int npts : number of points to use - def arc_function(self, r, rmin=0, rmax=None, use_points=False, base_cut=None): - rad, br, ct, center = self.fit_arc( - rmin=rmin, rmax=rmax, use_points=use_points) - #print('fitted radius, center, base_radius, cos(theta):',rad,center,br,ct) - return self.arc(r, center, rad) + return: + :tuple : (x,y,z) coordinates of points on the ellipsoid as ndarrays + """ + phi = np.arccos(2 * np.random.rand(npts) - 1) + theta = 2 * np.pi * np.random.rand(npts) + s = np.array([np.sin(phi) * np.cos(theta), np.sin(phi) + * np.sin(theta), np.cos(phi)]) + r = (np.array(parmsc['T'])@s).T + np.array(parmsc['v']).T + return r[:, 0], r[:, 1], r[:, 2] - @property - def contact_angle(self): - rad, br, ct, center = self.fit_arc() - return np.arccos(ct)*180./np.pi + @staticmethod + def _fit_circle(hr, hh, nonlinear=True): + """ fit an arc through the profile h(r) sampled by the class + :param list hr : list of arrays with the radial coordinates + :param list hh : list of arrays with the elevations + :param bool nonlinear : use the more accurate minimization of the rmsd instead of the algebraic distance + + return: + :list : a list with the tuple (radius, base radius, cos(theta), center, rmsd) + for each bin. If only one bin is present, return just the tuple. + """ + parms = [] + for i in np.arange(len(hr)): + r = hr[i] + h = hh[i] + if len(r) == 0: + parms.append(None) + break + M = np.vstack((r, h, np.ones(r.shape))).T + b = r**2 + h**2 + sol = np.linalg.lstsq(M, b, rcond=None)[0] + rc, hc = sol[:2]/2. + rad = np.sqrt(sol[2]+rc**2+hc**2) + pos = np.vstack([r, h]).T + if nonlinear: + res = minimize(ContactAngle.rmsd_circle, x0=[rad, rc, hc], args=(pos), + method='nelder-mead', options={'xatol': 1e-8, 'disp': False}) + rad, rc, hc = res.x + base_radius = np.sqrt(rad**2-hc**2) + rmsdval = res.fun + else: + rmsdval = ContactAngle.rmsd_circle([rad, rc, hc], pos) - @property - def base_radius(self): - rad, br, ct, center = self.fit_arc() - return br + base_radius = np.sqrt(rad**2-hc**2) + costheta = -hc/rad + theta = np.arccos(costheta) * 180./np.pi + if theta < 0: + theta += 180. - @property - def radius(self): - rad, br, ct, center = self.fit_arc() - return rad + parms.append((rad, base_radius, theta, [rc, hc], rmsdval)) + if len(parms) == 1: + return parms[0] + else: + return parms - @property - def center(self): - rad, br, ct, center = self.fit_arc() - return center - - def droplet_size(self, direction): - delta = np.max(self.inter.atoms.positions[:, direction]) - delta = delta - np.min(self.inter.atoms.positions[:, direction]) - return np.abs(delta) - - def remove_COM(self, removeCOM, droplet, inter, alpha, box): - if removeCOM is not None: - while(self.droplet_size(removeCOM) > box[removeCOM] - 4 * alpha): - pos = droplet.universe.atoms.positions.copy() - pos[:, removeCOM] += alpha - droplet.universe.atoms.positions = pos.copy() - pos = droplet.universe.atoms.pack_into_box() - droplet.universe.atoms.positions = pos.copy() - com = self.inter.atoms.center_of_mass()[removeCOM] - pos = droplet.universe.atoms.positions.copy() - pos[:, removeCOM] -= com - droplet.universe.atoms.positions = pos.copy() - - pos[:, removeCOM] += box[removeCOM]/2 - droplet.universe.atoms.positions = pos.copy() - - droplet.universe.atoms.positions = pos.copy() - pos = droplet.universe.atoms.pack_into_box() - droplet.universe.atoms.positions = pos.copy() - return droplet.universe.atoms.positions - - -class ContactAngleGitim(ContactAngle): - """ ContactAngle class implementation with GITIM - """ - - def sample(self, inter, bins=100, cut=3.5, alpha=2.5, pdbOut=None, binning='theta', periodic=None, removeCOM=None, base_cut=0.0): - """ compute the height profile z(r) of a droplet - - :param int bins : number of slices along the z axis (across the whole box z-edge) - :param float cut : cut-off for the clustering algorithm that separates liquid from vapor - :param float alpha : probe sphere radius for GITIM - :param str pdbOut : optional pdb file where to store trajectory with tagged surface atoms - - NOTES: - 1) not tested with the molecular option of GITIM - 2) bins in x,y and z directions are the same, this might be a problem for boxes with large aspect - ratios - 3) make removeCOM=0 to remove movement of droplet along x axis + @staticmethod + def _contact_angles_from_ellipse(p, off=0.0): + """ compute the contact angle from the parameters of the polynomial representation + + :parms array_like p : a sequence (list, tuple, ndarray) of parameters of the ellipse equation + in general form: + p[0] x^2 + p[1] x y + p[2] y^2 + p[3] x + p[4] y + p[5] = 0 + :param float off : elevation from the substrate surface, where the + contact angle should be evaluated. Default; 0.0, corresponding + to the atomic center of the highest atom of the substrate + return: + :tuple : a tuple (theta1,theta2) with the left and right (internal) + contact angles + + We require the solution of the line passing through the point to have two coincident + solutions for the system ax2 + bxy + cy2 + dx + ey +f = 0 and y-y0 = m (x-x0) + This yields m = 4*a*c*x0*y0 + 2*a*e*x0 - b**2*x0*y0 - b*d*x0 - b*e*y0 - 2*b*f + 2*c*d*y0 + d*e + The point x0,y0 is the one that solves the ellipse equation when requesting y0=off, so + x0 =[-b*y_0 - d +/- sqrt(-4*a*c*y_0**2 - 4*a*e*y_0 - 4*a*f + b**2*y_0**2 + 2*b*d*y_0 + d**2)]/(2*a) """ - droplet, substrate = self.droplet, self.substrate - box = droplet.universe.dimensions[:3] - try: - _ = self.r_ - except: - self.r_, self.z_, self.theta_, self.R_ = [], [], [], [] - self.nframes += 1 + a, b, c, d, e, f = p + if np.isclose(4*a*c, b**2, atol=1e-3): + rad = np.sqrt((d**2+e**2)/(4*a**2)-f/a) + hc = -e/(2*a) - off + base_radius = np.sqrt(rad**2-hc**2) + costheta = -hc/rad + theta = np.arccos(costheta) * 180./np.pi + if theta < 0: + theta += 180. + return (theta, theta) + y_0 = off + x_01 = (-b*y_0 - d + np.sqrt(-4*a*c*y_0**2 - 4*a*e*y_0 - + 4*a*f + b**2*y_0**2 + 2*b*d*y_0 + d**2))/(2*a) + x_02 = (-b*y_0 - d - np.sqrt(-4*a*c*y_0**2 - 4*a*e*y_0 - + 4*a*f + b**2*y_0**2 + 2*b*d*y_0 + d**2))/(2*a) + m1 = (4*a*c*x_01*y_0 + 2*a*e*x_01 - b**2*x_01*y_0 - + b*d*x_01 - b*e*y_0 - 2*b*f + 2*c*d*y_0 + d*e) + m1 = m1/(4*a*c*x_01**2 - b**2*x_01**2 - 2 * + b*e*x_01 + 4*c*d*x_01 + 4*c*f - e**2) + m2 = (4*a*c*x_02*y_0 + 2*a*e*x_02 - b**2*x_02*y_0 - + b*d*x_02 - b*e*y_0 - 2*b*f + 2*c*d*y_0 + d*e) + m2 = m2/(4*a*c*x_02**2 - b**2*x_02**2 - 2 * + b*e*x_02 + 4*c*d*x_02 + 4*c*f - e**2) + # depending on the sign of a at the denominator of x_01 and x_02, they can have a different order + # along the axis: let's keep them ordered, so that the first is the left one and the second the right one. + if x_02 < x_01: + x_01, x_02 = x_02, x_01 + m1, m2 = m2, m1 - self.inter = inter + theta1 = np.arctan(m1)*180/np.pi + theta2 = np.arctan(m2)*180/np.pi + # we compute the internal angle (assuming the droplet is in the upper half space), and need to take care of this + # theta1 is at the left edge of the droplet, theta2 at the right + if theta1 < 0.0: + theta1 += 180. + if theta2 < 0.0: + theta2 = -theta2 + else: + theta2 = 180.-theta2 - self.remove_COM(removeCOM, droplet, inter, alpha, box) - self.maxr = np.nanmax(droplet.universe.dimensions[:2]) - self.maxz = droplet.universe.dimensions[2] + return (theta1, theta2) - if pdbOut is not None: - self.inter.writepdb(pdbOut, multiframe=True) + @staticmethod + def _fit_ellipse_fitzgibbon(x, y): + D = np.vstack([x**2, x*y, y**2, x, y, np.ones(len(x))]).T + S = D.T @ D + C = np.zeros((6, 6), dtype=float) + C[0, 2], C[2, 0], C[1, 1] = 2, 2, -1 + eigval, eigvec = np.linalg.eig(np.linalg.solve(S, C)) + sort = np.argsort(eigval) + eigval, eigvec = eigval[sort], eigvec[sort] + lam = np.nonzero(np.logical_and( + eigval > 0, np.isfinite(eigval)))[0][-1] + return eigvec[lam] - if periodic is None: - r, theta, z = self.spherical_coordinates() - else: - r, theta, z = self.cylindrical_coordinates(periodic) + @staticmethod + def _fit_ellipse_flusser(x, y): + D1 = np.vstack([x**2, x*y, y**2]).T + D2 = np.vstack([x, y, np.ones(len(x))]).T + + S1 = D1.T @ D1 + S2 = D1.T @ D2 + S3 = D2.T @ D2 + T = -np.linalg.inv(S3) @ S2.T + M = S1 + S2 @ T + C = np.array(((0, 0, 2), (0, -1, 0), (2, 0, 0)), dtype=float) + M = np.linalg.solve(C, M) + eigval, eigvec = np.linalg.eig(M) + sort = np.argsort(eigval)[::-1] + eigval, eigvec = eigval[sort], eigvec[sort] + con = 4 * eigvec[0] * eigvec[2] - eigvec[1]**2 > 0 + ak = eigvec[:, np.nonzero(con)[0]] + return np.concatenate((ak, T @ ak)).ravel() - if binning == 'theta': - self.sample_theta_R(bins, theta, r, z, base_cut) - elif binning == 'z': - self.sample_z_r(bins, z, r, base_cut) + @staticmethod + def _fit_ellipsoid_lstsq_cox(x, y, z): + """ + Return the polynomial coefficients that perform a least square + fit to a set of points, using a modified version of: + Turner, D. A., I. J. Anderson, J. C. Mason, and M. G. Cox. + "An algorithm for fitting an ellipsoid to data." + National Physical Laboratory, UK (1999). + + :param array_like x : list or ndarray array with coordinate x + :param array_like y : list or ndarray array with coordinate y + :param array_like z : list or ndarray array with coordinate z + + return: + :ndarray : a (10,)-ndarray with the polynomial coefficients + """ + D = np.vstack([x**2+y**2-2*z**2, x**2 - 2*y**2 + z**2, + 4*x*y, 2*x*z, 2*y*z, x, y, z, np.ones(len(x))]).T + E = np.vstack([x**2+y**2+z**2]).T + DTD, DTE = D.T@D, D.T@E + p2 = scipy.linalg.solve(DTD, DTE).flatten() + u, v, m, n, p, q, r, s, t = p2 + a, b = (1-u-v), (1-u+2*v) + c = (3.-a-b) + # a-c = a- 3 + a + b = 2a + b - 3 = 2-2u-2v +1 -u+2v -3 = -3 u + # a-b = 1-u-v -1 +u -2 v = -3 v + d, e, f = -2*m, -n, -p + g, h, k = -q, -r, -s + l = -t + p = a, b, c, f, e, d, g, h, k, l + if any(np.array([4*a*b-d*d, 4*a*c-e*e, 4*b*c-f*f]) < 0): + print('non positive definite') + return np.array(p) + + @staticmethod + def _fit_ellipse(x, y, nonlinear=True, off=0.0, points_density=25, verbose=False): + """ fit an ellipse through the points (x,y) + + :param list x : list of arrays with coordinates x + :param list y : list of arrays with coordinates y + :param bool nonlinear : use the more accurate minimization of the rmsd instead of the algebraic distance + :param float off : elevation from the substrate surface, where the + contact angle should be evaluated. Default; 0.0, corresponding + to the atomic center of the highest atom of the substrate + :param int points_density: number of points per Angstrom on the ellipse that are used to compute the rmsd + + return: + + :list : a list with a tuple (parms,parmsc,theta, rmsd) for each of the bins + If only one bin is present, return just the tuple + parms : array of parameters of the ellipse equation in general form: + ellipse: a[0] x^2 + a[1] x y + a[2] y^2 + a[3] x + a[4] y = 0 + parmsc : dictionary with parameters of the ellipse canoncial form: (x0,y0,a,b,phi,e) + with a,b the major and minor semiaxes, x0,y0 the center and theta the angle + between x axis and major axis + theta : [left angle, right angle] + rmsd : the rmsd to the best fit (linear or nonlinear) ellipse or ellipsoid + + Uses a slighly modified version of Fitzgibbon's least square algorithm from + Stable implementation from Halır, Radim, and Jan Flusser, + "Numerically stable direct least squares fitting of ellipses." + Proc. 6th International Conference in Central Europe on Computer + Graphics and Visualization. WSCG. Vol. 98. Citeseer, 1998. + + python code for ellipse fitting from https://scipython.com/blog/direct-linear-least-squares-fitting-of-an-ellipse/ + """ + retvals = [] + + for i in np.arange(len(x)): + XX = x[i] + YY = y[i] + if len(XX) == 0: + retvals.append(None) + break + idx1 = np.isnan(XX) + idx2 = np.isnan(YY) + idx = np.logical_and(~idx1, ~idx2) # nan removed + X, Y = XX[idx], YY[idx] + p = ContactAngle._fit_ellipse_flusser(X, Y) + pos = np.vstack([X, Y]).T + # we aim at accuracy of about 0.04 Angstrom with points_density = 25 + N = int(points_density * np.max(np.sqrt(X**2+Y**2))) + if len(p) == 0: + # no positive eigenval, we start over from the (slightly perturbed) circle + rad, rb, costheta, [x0, y0], rmsd = ContactAngle._fit_circle( + [X], [Y], nonlinear=nonlinear) + # set a=c=1 in Eqs 27-30 https://mathworld.wolfram.com/Circle.html + p = [1.0001, 0.0, 0.9999, -2*x0, -2*y0, x0**2+y0**2-rad**2] + if nonlinear: + rmsdf = ContactAngle.rmsd_ellipse_penalty + res = minimize(rmsdf, x0=p, args=(pos, N, False), + method='nelder-mead', options={'xatol': 1e-4, 'disp': verbose}) + p, rmsdval = res.x, res.fun + else: + rmsdval = ContactAngle.rmsd_ellipse(p, pos, N) + theta1, theta2 = ContactAngle._contact_angles_from_ellipse( + p, off=off) + # canonical form (careful: we overwrite internal variables a,b) + cp = ContactAngle._ellipse_general_to_canonical( + p, check_coeffs=False) + retvals.append((p, cp, [theta1, theta2], rmsdval)) + if len(retvals) == 1: + return retvals[0] else: - raise (ValueError("Wrong binning type")) + return retvals + + @staticmethod + def _fit_ellipsoid(x, y, z, nonlinear=True, off=0.0, points_density=25): + """ fit an ellipsoid through the points (x,y,z) + + :param list x : list of arrays with coordinates x + :param list y : list of arrays with coordinates y + :param list z : list of arrays with coordinates z, + :param bool nonlinear : use the more accurate minimization of the rmsd instead of the algebraic distance + :param float off : elevation from the substrate surface, where the + contact angle should be evaluated. Default; 0.0, corresponding + to the atomic center of the highest atom of the substrate + :param int points_density: number of points per Angstrom on the ellipse that are used to compute the rmsd + + return: + + :list : a list with a tuple (parms,parmsc,theta, rmsd) for each of the bins + If only one bin is present, return just the tuple + parms : array of parameters of the ellipsoid equation in general form: + a[0] x^2 + a[1] y^2 + a[2] z^2 + a[3] yz + a[4] xz + + a[5] xy + a[6] x + a[7] y + a[8] z + a[9] = 0, otherwise. + parmsc : dictionary with parameters of the ellipsoid affine form (T,v), + such that the ellipsoid points are + r = T s + v, + if s are points from the unit sphere centered in the origin. + theta : the contact angle as a function of the azimuthal angle phi from phi=0, aligned + with (-1,0,0) to 2 pi. + rmsd : the rmsd to the best fit (linear or nonlinear) ellipsoid + + For the least square approach see _fit_ellipsoid_lstsq_cox - def _compute_coords(self, z_, r_): - # r_,distance from center;R_, projection on xy plane can be considered as x in 2d - R_ = np.sqrt(r_*r_ + z_*z_) - theta_ = np.arccos(r_/R_) - self.theta_ += list(theta_) - self.r_ += list(r_) - self.z_ += list(z_) - #self.R_ += list(R_) + """ + retvals = [] + for i in np.arange(len(x)): + XX, YY, ZZ = x[i], y[i], z[i] + if len(XX) == 0: + retvals.append(None) + break + idx1, idx2, idx3 = np.isnan(XX), np.isnan(YY), np.isnan(ZZ) + idx = np.logical_and(~idx1, ~idx2) + idx = np.logical_and(idx, ~idx3) # nan removed + X, Y, Z = XX[idx], YY[idx], ZZ[idx] + p = ContactAngle._fit_ellipsoid_lstsq_cox(X, Y, Z) + pos = np.vstack([X, Y, Z]).T + # we aim at accuracy of about 0.04 Angstrom with points_density = 25 + N = int(points_density * np.max(np.sqrt(X**2+Y**2))) + if nonlinear: + rmsdf = ContactAngle.rmsd_ellipsoid_penalty + start = rmsdf(p, pos, N, False) + # we aim at improving by 10% the least square fit + res = minimize(rmsdf, x0=p, args=(pos, N, False), + method='nelder-mead', options={'xatol': start/10, 'disp': False}) + p, rmsdval = res.x, res.fun + else: + rmsdval = ContactAngle.rmsd_ellipsoid(p, pos, N) - r, z, theta = np.asarray(self.r_), np.asarray( - self.z_), np.asarray(self.theta_) + theta = np.array([]) + cp = ContactAngle._ellipsoid_general_to_affine( + p, check_coeffs=False) + retvals.append((p, cp, theta, rmsdval)) - return r, theta, z + if len(retvals) == 1: + return retvals[0] + else: + return retvals - def cylindrical_coordinates(self, periodic): - pos = self.remove_com(self.inter.atoms) - dirs = np.array([0, 1]) - dd = dirs[dirs != periodic][0] - r_ = pos[:, dd] # change to this if you face error "np.abs(pos[:,dd])" - z_ = pos[:, 2] - # R_=np.sqrt(R_*R_+z_*z_) - # print(z_,r_) - # print(r_.shape,z_.shape) - return self._compute_coords(z_, r_) - - def spherical_coordinates(self): - pos = self.remove_com(self.inter.atoms) - #r_ = np.linalg.norm(pos[:,0:2],axis=1) - # this r_ is the distance from the center - r_ = np.linalg.norm(pos[:, 0:2], axis=1) - z_ = pos[:, 2] - return self._compute_coords(z_, r_) - - def sample_theta_R(self, bins, theta, r, z, base_cut): - R = np.sqrt(r*r+z*z) - n_h, t_edges = np.histogram(theta, bins=bins, range=( - 0, np.pi), weights=None, density=False) - R_h, t_edges = np.histogram(theta, bins=bins, range=( - 0, np.pi), weights=R, density=False) - #cond = np.where(n_h>0) - zz, rr = R_h * np.sin(t_edges[:-1]) / \ - (1.*n_h), R_h * np.cos(t_edges[:-1])/(1.*n_h) - zzmin = np.nanmin(zz) - ercut = zzmin+base_cut - self.z, self.r = zz[zz > ercut], rr[zz > ercut] - # print("sample",z,r) - - def sample_z_r(self, bins, z, r, base_cut): - #print("I love maxr",self.maxr) - n_h, r_edges = np.histogram(r, bins=bins, range=( - 0, self.maxr), weights=None, density=False) - z_h, r_edges = np.histogram(r, bins=bins, range=( - 0, self.maxr), weights=z, density=False) - zz, rr = z_h/(1.*n_h), r_edges[:-1] - zzmin = np.nanmin(zz) - ercut = zzmin+base_cut - self.z, self.r = zz[zz > ercut], rr[zz > ercut] - #self.z , self.r = z_h/(1.*n_h), r_edges[:-1] - - def fit_arc(self, rmin=0, rmax=None, use_points=False): - """ fit an arc through the profile z(r) sampled by the class - - :param float rmin : minimum radius used for the fit (default: 0) - :param float rmax : maximum radius used for the fit (default: maximum from dataset) - :param bool use_points : False: use the mean surface atoms positions - True : use the surface atoms positions + @staticmethod + def _ellipse_canonical_to_general(coeffs): + """ Convert canonical coefficients (x0,y0,A,B,phi) to general ones (a,b,c,d,e,f) """ + x0, y0, A, B, phi = coeffs + a = A*A * np.sin(phi)**2 + B*B*np.cos(phi)**2 + b = 2*(A*A-B*B)*np.sin(phi)*np.cos(phi) + c = A*A * np.cos(phi)**2 + B*B*np.sin(phi)**2 + d = -2*a*x0 - b*y0 + e = -b*x0 - 2*c*y0 + f = a*x0**2 + b * x0*y0 + c*y0**2 - a*a*b*b + return [a, b, c, d, e, f] - try: - if use_points: - z, r = self.z_, self.r_ - else: - z, r = self.z, self.r + @staticmethod + def _ellipse_general_to_canonical(coeffs, check_coeffs=True): + """ Convert general coefficients (a,b,c,d,e,f) to canonical ones. - except AttributeError: - raise RuntimeError( - 'something is wrong: no surface atoms or not gibbs surface present') + :param list coeffs : general coefficients + :param bool check_coeffs : raise an error if the coefficients do not represent + an ellipse. In some cases this checks needs to be + turned off, e.g. for a COBYLA minimization, as + during the minimization the constraints might + be violated. - return self.fit_arc_(z, r, rmin=rmin, rmax=rmax) + return: + : dict : a dictionary with the canonical coefficients x0,y0,a,b,phi,e - def fit_points_ellipse(self, rmax=None, rmin=None, bins=None, yy=None, use_points=False): - if rmax is None: - rmax, rmin, bins = 80, -80, 3200 - yy = 0 - try: - if use_points: - z, r = self.z_, self.r_ + from https://scipython.com/blog/direct-linear-least-squares-fitting-of-an-ellipse/ + + Convert the cartesian conic coefficients, (a, b, c, d, e, f), to the + ellipse parameters, where F(x, y) = ax^2 + bxy + cy^2 + dx + ey + f = 0. + The returned parameters are x0, y0, ap, bp, e, phi, where (x0, y0) is the + ellipse centre; (ap, bp) are the semi-major and semi-minor axes, + respectively; e is the eccentricity; and phi is the rotation of the semi- + major axis from the x-axis. + + """ + # We use the formulas from https://mathworld.wolfram.com/Ellipse.html + # which assumes a cartesian form ax^2 + 2bxy + cy^2 + 2dx + 2fy + g = 0. + # Therefore, rename and scale b, d and f appropriately. + a, b, c, d, f, g = coeffs + b /= 2 + d /= 2 + f /= 2 + den = b**2 - a*c + if check_coeffs and den > 0: + raise ValueError('coeffs do not represent an ellipse: b^2 - 4ac must' + ' be negative!'+str(den)) + # The location of the ellipse centre. + x0, y0 = (c*d - b*f) / den, (a*f - b*d) / den + num = 2 * (a*f**2 + c*d**2 + g*b**2 - 2*b*d*f - a*c*g) + fac = np.sqrt((a - c)**2 + 4*b**2) + # The semi-major and semi-minor axis lengths (these are not sorted). + if np.isclose(fac, 0.0): # a=c and b=0 + ap, bp = a, a + else: + argument = num / den / (fac - a - c), num / den / (-fac - a - c) + if den > 0: # it is not an ellipse, so not all semi axes are defined. + # for the purpose of the minimizer, we still need to provide + # one solution + ap, bp = np.sqrt([np.max(argument)]*2) else: - z, r = self.z, self.r - except AttributeError: - raise RuntimeError( - 'something is wrong: no surface atoms or not gibbs surface present') - cofX, th1, th2 = self.fit_ellipse_(r, z, yy) - - print(cofX, th1, th2) - return self.arc_ellipse(cofX, rmax=rmax, rmin=rmin, bins=bins) - - def sample_theta_R(self, bins, th_left, th_right, RR_left, RR_right): - n_h_left, t_edges_left = np.histogram( - th_left, bins=bins, range=(0, np.pi), weights=None, density=False) - R_h_left, t_edges_left = np.histogram(th_left, bins=bins, range=( - 0, np.pi), weights=RR_left, density=False) - n_h_right, t_edges_right = np.histogram( - th_right, bins=bins, range=(0, np.pi), weights=None, density=False) - R_h_right, t_edges_right = np.histogram(th_right, bins=bins, range=( - 0, np.pi), weights=RR_right, density=False) - self.left_z, self.left_r = R_h_left * np.sin(t_edges_left[:-1])/(1.*n_h_left),\ - R_h_left * np.cos(t_edges_left[:-1])/(1.*n_h_left) - self.right_z, self.right_r = R_h_right * np.sin(t_edges_right[:-1])/(1.*n_h_right),\ - R_h_right * np.cos(t_edges_right[:-1])/(1.*n_h_right) - - -class ContactAngleGibbs(ContactAngle): - """ ContactAngle class implementation using an estimate of the - Gibbs dividing surface - """ - - def sigmoid(self, r, r0, A, w): - return A*(1.+np.tanh((r0-r)/w))/2. - - def sample(self, bins=100, params=None, binning='theta', periodic=None, base_cut=0.0): - """ compute the height profile z(r) of a droplet - :param int bins : number of slices along the z axis (across the whole box z-edge) + ap, bp = np.sqrt(argument) + # Sort the semi-major and semi-minor axis lengths but keep track of + # the original relative magnitudes of width and height. + width_gt_height = True + if ap < bp: + width_gt_height, ap, bp = False, bp, ap + # The eccentricity. + r = (bp/ap)**2 + if r > 1: + r = 1./r + e = np.sqrt(1. - r) + # The angle of anticlockwise rotation of the major-axis from x-axis. + if np.isclose(b, 0) or np.isclose(ap, bp): # handle also the circle + phi = 0.0 if a <= c else np.pi/2. + else: + phi = np.arctan((2.*b) / (a - c)) / 2. + if a > c: + phi += np.pi/2. + if not width_gt_height: + # Ensure that phi is the angle to rotate to the semi-major axis. + phi += np.pi/2. + # NOTE: don't change the units of phi to deg., other parts of the code + # depend on it being in rad. + phi = phi % np.pi + return {'x0': x0, 'y0': y0, 'a': ap, 'b': bp, 'phi': phi, 'e': e} + + @staticmethod + def _ellipsoid_base(coeffs, off, npts=1000): + # sets z=off and expresses the polynomial coefficients + # of an ellipse in terms of the intersection with the ellipsoid + # then calculates the points in the Cartesian plane. + a, b, c, f, g, h, p, q, r, d = coeffs + coeffs_ellipse = a, 2*h, b, 2*g*off+p, 2*f*off+q, r*off+d + pc = ContactAngle._ellipse_general_to_canonical(coeffs_ellipse) + x, y = ContactAngle.ellipse(pc, npts=npts, tmin=0, tmax=2.*np.pi) + return x, y + + @staticmethod + def _ellipsoid_normal(coeffs, x, y, z): + # first we need to recover the values of (theta,phi) corresponding to the + # points x,y,z belonging to the ellipsoid defined by the polynomial with + # coefficients coeffs + T, v = ContactAngle._ellipsoid_general_to_affine(coeffs).values() + pos = np.vstack([x, y, z]).T + # points on the unit sphere + s = np.dot(np.linalg.inv(T), (pos - v).T).T + s = s.T / np.linalg.norm(s, axis=1) # tidy up roundoff errors + sx, sy, sz = s[0], s[1], s[2] + theta = np.arccos(sz) # radius is one + phi = np.sign(sy) * sx / np.sqrt(1-sz**2) + dsdtheta = np.vstack( + [np.cos(theta)*np.cos(phi), np.cos(theta)*np.sin(phi), -np.sin(theta)]).T + dsdphi = np.vstack([-np.sin(theta)*np.sin(phi), + np.sin(theta)*np.cos(phi), 0*phi]).T + normal = np.cross(dsdtheta, dsdphi) + # we need now to check if it's pointing outwards or inwards. + # It's enough to inspect one point: since v is the center of + # the affine transformation, x-v is the displacement from it, + # and we request that (x-v).n > 0 for n to be pointing out. + if (pos[0]-v).T@normal[0] < 0: + normal = -normal + # tidy up roundoff errors + return (normal.T/np.linalg.norm(normal, axis=1)).T + @staticmethod + def _ellipsoid_contact_line_and_angles(coeffs, off, npts=1000): + """ Computes the contact line and the contact angles along the + line. + + :param array_like coeffs : the coefficients of the polynomial + defining the ellipsoid, see, e.g., + _ellipsoid_general_to_affine + :param float off : elevation above the topmost atoms + of the substrate, used to compute + the elliptical intersection that is + the contact line + :param int npts : number of points to sample + + return: + :tuple : a tuple containing the contact line + on the (x,y) plane as a (npts,2)-ndarray + and the contact angles as a (npts,)-ndarray. """ - from scipy.optimize import curve_fit + bx, by = ContactAngle._ellipsoid_base(coeffs, off, npts=npts) + bz = np.zeros(bx.shape)+off + n = ContactAngle._ellipsoid_normal(coeffs, bx, by, bz) + theta = np.arccos(n[:, 2]) + return (np.vstack([bx, by]).T, theta) - droplet, substrate = self.droplet, self.substrate - self.nframes += 1 + @staticmethod + def _ellipsoid_general_to_affine(coeffs, check_coeffs=True): + """ Convert general coefficients (a,b,c,f,g,h,p,q,r,d) of the polynomial + a xx + b yy +c zz + 2f yz + 2g xz + 2h xy + px + qy +rz +d to the affine + transformation ones that map the unitary sphere centered in the origin + to a point on the ellipsoid. + + :param list coeffs : general coefficients + :param bool check_coeffs : raise an error if the coefficients do not represent + an ellipse. In some cases this checks needs to be + turned off, e.g. for a COBYLA minimization, as + during the minimization the constraints might + be violated. + + return: + : dict : a dictionary with the affine matrix and vector, T and v + If a point on the unit sphere centered in the origin is s, + then the correspoinding point r=(x,y,z) on the ellipsoid is + defined by the affine transofmation r = T s + v = A^(-1/2) s + v + + Here, A is the matrix of the quadratic form defining the ellipsoid, + + r^t A r + beta^t r + gamma = 0. + + The connection to the general coefficients is given by + + | a h g | + A= | h b f | , beta^t = (p, q, r) , gamma = d + | g f c | + + In addition, the quadratic form can be recast to (r-v)^t A' (r-v) = k. + This form is not unique, but this is helpful to find the correspondence + to the polynmomial coefficients, which are also defined up to an immaterial + multiplicative constant. Note that setting k=1 imposes a constaint on gamma, + but gamma=d is given, so we need to keep the general form and determine k. + After expansion and comparison, one finds + + A' = A , v = -1/2 A^-1 beta , k = (1/4) beta^t A^-1 beta - gamma, + + Now that k is known, we can recast the equation in the standard form + (r-v)^t M (r-v) = 1 by setting M = A'/k = A/k. Note that v has to be determined + using A and not A/k. The standard form can be used to generate points r on the + ellipse starting from the unit sphere s centered in zero as the affine transformation + + r = v + T s, with T = A^{-1/2} + + This function eventually returns the elements of the affine transofmation T and v. + + >>> # This tests that we can recover the original parameters from a fit + >>> p1 = np.array([2, 3, 4., 0.4, 0.5, -0.6, 0.3, 4.2, 0.1, -10]) + >>> cp = ContactAngle._ellipsoid_general_to_affine(p1,check_coeffs=True) + >>> x,y,z = ContactAngle.ellipsoid(cp) + >>> p2 = ContactAngle._fit_ellipsoid_lstsq_cox(x,y,z) + >>> # NOTE: for Cox, a+b+c = 3, so, to compare the two sets we need to normalize + >>> p2 *= np.sum(p1[:3])/3. + >>> print(np.all(np.isclose(p1,p2))) + True + """ + + a, b, c, f, g, h, p, q, r, d = coeffs + if check_coeffs: + conds = [a*b - h ** 2 > 0, a*c - g**2 > 0, b*c - f**2 > 0] + if not all(conds): + raise ValueError('coeffs do not represent an ellipsoid!') + A = np.array([[a, h, g], [h, b, f], [g, f, c]]) + beta = np.array([p, q, r]) + norm = -d + beta.T@np.linalg.inv(A)@beta/4 + A /= norm + U, Sigma, VT = np.linalg.svd(A) + invsqrt_A = np.linalg.inv(U@np.diag(np.sqrt(Sigma))@VT) + v = -0.5 * np.linalg.inv(A*norm)@beta + return {'T': invsqrt_A.tolist(), 'v': v.tolist()} + + def _compute_coords(self, r, h, symm): + # r is the distance of the point projected on the + # substrate from the axis/center of the + # cylindrical/sperical droplet + # + # h is the elevation from the substrate surface + # + # symm is the remaining coordinate (the angle phi in + # spherical coordinates, the position along the cylinder + # axis in cylindrical ones. + + R = np.sqrt(r**2 + h**2) + theta = np.arccos(r/R) + return np.asarray(r), np.asarray(theta), np.asarray(h), np.asarray(symm) - pos = self.remove_com(droplet) - if periodic is None: - r = np.linalg.norm(pos[:, 0:2], axis=1) - z = pos[:, 2] - R = np.sqrt(r*r + z*z) + def cartesian_coordinates(self, group): + pos = self._shifted_pos(group) + if self.periodic is None: + return pos[:, 0], pos[:, 1], pos[:, 2] else: dirs = np.array([0, 1]) - dd = dirs[dirs != periodic][0] - r = np.abs(pos[:, dd]) - z = pos[:, 2] - R = np.sqrt(r*r+z*z) + dd = dirs[dirs != self.periodic][0] + od = list(set((1, 2, 3))-set((2, dd)))[0] + return pos[:, dd], pos[:, od], pos[:, 2] - theta = np.arccos(r/R) - maxr = np.nanmax(droplet.universe.dimensions[:2]) - maxR = maxr - maxz = droplet.universe.dimensions[2] - if binning == 'z': - # histogram for this frame - H_, zedge, redge = np.histogram2d( - z, r, bins, [[0, maxz], [0, maxr]]) - zedge = (zedge[:-1]+zedge[1:])/2. - redge = (redge[:-1]+redge[1:])/2. - # print(H_[H_>0],zedge,redge) - elif binning == 'theta': - H_, thetaedge, Redge = np.histogram2d( - theta, R, bins, [[0, np.pi/2.], [30, maxR]]) - thetaedge = (thetaedge[:-1]+thetaedge[1:])/2. - Redge = (Redge[:-1]+Redge[1:])/2. + def cylindrical_coordinates(self, group): + pos = self._shifted_pos(group) + dirs = np.array([0, 1]) + dd = dirs[dirs != self.periodic][0] + r = pos[:, dd] + h = pos[:, 2] + od = list(set((1, 2, 3))-set((2, dd)))[0] + z = pos[:, od] + return self._compute_coords(r, h, z) + + def spherical_coordinates(self, group): + pos = self._shifted_pos(group) + r = np.linalg.norm(pos[:, 0:2], axis=1) + h = pos[:, 2] + symm = np.arctan2(pos[:, 1], pos[:, 0]) + phi = np.arctan2(h, r)*180/np.pi + phi[phi < 0] += 360. + return self._compute_coords(r, h, phi) + + def sample_theta_R(self, theta, r, h): + """ + given a set of angles (theta), radii (r) and elevations (h), compute the mean profile h(r) + by taking the binned average values of h and r. + """ + R = np.sqrt(r*r+h*h) + n_h, t_edges = np.histogram( + theta, bins=self.bins, density=False, range=(0, np.pi), weights=None) + R_h, t_edges = np.histogram( + theta, bins=self.bins, density=False, range=(0, np.pi), weights=R) + cond = n_h > 0.0 + R_h = R_h[cond] + t_edges = t_edges[:-1][cond] + n_h = n_h[cond] + hh, rr = R_h * np.sin(t_edges)/n_h, R_h * np.cos(t_edges)/n_h + hhmin = np.min(hh) + ercut = hhmin+self.hcut + cond = hh > ercut + self.histo.x, self.histo.h = rr[cond], hh[cond] + + def _select_coords(self, use, bins): + # return lists of array with coordinates falling into the bins that partition the symmetric coordinate. + # The main idea is to use the instantaneous or stored cartesian coordinates (renamed so that z is always + # the surface normal and, in case of cylyndrical symmetry, y points along the cylinder axis. + # The cae of the histogram is different, because the relation, say (x,z) is not necessarily a function. + # In this case, we store the spherical/cylindrical coordinate values in the histogram, and we recover the + # averaged location of the surface a posteriori. + if self._sampled == False: + self._sampled, _ = True, self.sample() + try: + if self.store and use == 'stored': + x, h = self._x, self._h + if self.periodic is None: + y = self._y + elif not self.store and use == 'stored': + raise (ValueError( + "The choice use='stored' can only be used if store=True was passed at initialization ")) + elif use == 'frame': + x, h = self.x, self.h + if self.periodic is None: + y = self.y + elif use == 'histogram': + x, h = self.histo.x, self.histo.h + if self.periodic is None: + y = self.histo.y + else: + raise (ValueError( + "The parameter use can only take the values 'frame', 'histogram', or 'stored' ")) + + except AttributeError: + raise RuntimeError('No surface atoms or Gibbs surface present') + + # we always comply with the request of cutting all molecules below hcut from the minimum + if self.hcut > 0: + hmin = np.min(h) + cond = h > hmin+self.hcut + x, h = x[cond], h[cond] + if self.periodic is None: + y = y[cond] + if bins > 1: + # TODO: extend for histogram and stored + if use != 'frame': + raise (ValueError("bins>1 can be used only with use='frame'")) + hx, hy, hh = [], [], [] + symm = self.phi[cond] if self.periodic is None else self.y[cond] + limit = 2 * \ + np.pi if self.periodic is None else self.substrate.dimensions[self.periodic] + # set the right edges of the bins. This assumes that shifted_pos() has been called on the + # coordinates (hence, this function should not be public) + binvals = np.linspace(-limit/2., limit/2., bins+1)[1:-1] + # this will put in bin index 0 eveything close to or below 0.0 and in bin index + # nbins-1 everyhing close to or above limit + inds = np.digitize(symm, binvals) + for i in range(bins): + hx.append(x[inds == i]) + hh.append(h[inds == i]) + if self.periodic is None: + hy.append(y[inds == i]) else: - raise ValueError("binning can be only 'z' or 'theta'") + hx, hh = [x], [h] + if self.periodic is None: + hy = [y] + if self.periodic is None: + return hx, hy, hh + else: + return hx, hh - # cumulative histogram - try: - self.H - self.H += H_ + def fit_circle(self, use='frame', nonlinear=True, bins=1): + """ fit an arc through the profile h(r) sampled by the class - except: - self.H = H_ + :param str use : 'frame' : use the positions of the current frame only (default) + 'histogram': use the binned values sampled so far + 'stored' : use the stored surface atoms positions, if the option store=True was passed at initialization + :param int bins : the number of bins to use along the symmetry direction (cylinder axis, azimuthal angle) - # normalize by number of frames and metric factor + return: + a list including, for each of the bins: + : tuple : radius, base radius, cos(theta), center + """ + r, h = self._select_coords(use, bins=bins) + + return self._fit_circle(r, h, nonlinear=nonlinear) + + def fit_ellipse(self, use='frame', nonlinear=True, bins=1): + """ fit an ellipse through the points sampled by the class. See implementation details in _fit_ellipsoid() + + :param str use : 'frame' : use the positions of the current frame only (default) + 'histogram': use the binned values sampled so far + 'stored' : use the stored surface atoms positions, if the option store=True + was passed at initialization + :param int bins : the number of bins to use along the symmetry direction (cylinder axis, azimuthal angle) + + :return: + a list including, for each of the bins, a tuple with elements: + parms : parameters of the ellipse polynomial in general form: + a[0] x^2 + a[1] x y + a[2] y^2 + a[3] x + a[4] y = 0 + parmsc : dictionary of parameters in canoncial form: (a,b, x0,y0,phi, e) + with a,b the major and minor semiaxes, x0,y0 the center, phi the angle + (in rad) between x axis and major axis, and e the eccentricity. + theta : [left contact angle, right contact angle] + """ - # fit the histogram with a sigmoidal function, determine a Gibbs dividing distance for each z-slice - # we redo this every frame, it does not take much time, and we always have the up-to-date location - # of the Gibbs dividing surface with maximum statistics - if binning == 'z': - if periodic is None: - H_ = self.H / (self.nframes*redge) - else: - H_ = self.H / (self.nframes) # need to check + r, h = self._select_coords(use, bins=bins) + p, cp, theta, rmsd = self._fit_ellipse( + r, h, nonlinear=nonlinear, off=0.0) + self._polynomial_coefficients, self._canonical_form = p, cp + return p, cp, theta, rmsd + + def fit_ellipsoid(self, use='frame', nonlinear=True, bins=1): + """ fit an ellipsoid through the points sampled by the class. See implementation details in _fit_ellipsoid() + + :param str use : 'frame' : use the positions of the current frame only (default) + 'histogram': use the binned values sampled so far + 'stored' : use the stored surface atoms positions, if the option store=True + was passed at initialization + :param int bins : the number of bins to use along the symmetry direction (cylinder axis, azimuthal angle) + + :return: + :list : a list with a tuple (parms,parmsc,theta, rmsd) for each of the bins + If only one bin is present, return just the tuple + parms : array of parameters of the ellipsoid equation in general form: + a[0] x^2 + a[1] y^2 + a[2] z^2 + a[3] yz + a[4] xz + + a[5] xy + a[6] x + a[7] y + a[8] z + a[9] = 0, otherwise. + parmsc : dictionary with parameters of the ellipsoid affine form (T,v), + such that the ellipsoid points are + r = T s + v, + if s are points from the unit sphere centered in the origin. + theta : the contact angle as a function of the azimuthal angle phi from phi=0, aligned + with (-1,0,0) to 2 pi. + rmsd : the rmsd to the best fit (linear or nonlinear) ellipsoid + """ - elif binning == 'theta': - if periodic is None: - H_ = self.H / (self.nframes*Redge**2 * np.cos(thetaedge)) + x, y, z = self._select_coords(use, bins=bins) # TODO FIXME - else: - H_ = self.H / (self.nframes*Redge * np.cos(thetaedge)) # check - #print(self.H, H_, self.nframes) - - parms, zvals, rvals = [], [], [] - for i, h in enumerate(H_): - cond = h > 1e-6 - if params is None: - p0 = (40., 10.0, 0.8) - else: - p0 = params - try: - if binning == 'z': - popt, pcov = curve_fit( - self.sigmoid, redge[cond], h[cond], p0=p0) - #popt, pcov = curve_fit(self.exp_sigmoid, redge[cond], h[cond], p0=p0) - - if binning == 'theta': - popt, pcov = curve_fit( - self.sigmoid, Redge[cond], h[cond], p0=p0) - - except TypeError: - pass # not enough points in slice i - except ValueError: - # not enough points (another version of numpy throws this error? check) - pass - except RuntimeError: - pass # fit failed - try: - if np.isfinite(pcov[0][0]) and popt[0] != p0[0] and popt[0] > 0 and popt[0] < maxR: - #print(pcov[0][0], popt[0]) - # if np.isfinite(pcov[0][0]) and popt[0]>0 and popt[0] flcut], rvals[zvals > flcut] + p, cp, theta, rmsd = self._fit_ellipsoid( + x, y, z, nonlinear=nonlinear, off=0.0) + self._polynomial_coefficients, self._canonical_form = p, cp + return p, cp, theta, rmsd