Merge pull request #67 from WildJimmy/main

Added Jupyter Notebook for Configurations 3 and 4, as well as Position Finding .py File

position_guess.py

@@ -0,0 +1,129 @@
+'''
+This is the "brute force" python file I was using to find/test
+different initial positions, which I've added here for completeness.
+It's sloppy but (I think) it was able to converge on most of them given the
+proper filters
+'''
+
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import itertools
+import scipy.optimize
+
+
+ke2 = 197 / 137  # eV-nm   Coulomb force charge
+alpha = 1.09e3  # eV      parameter of model
+rho = 0.0321  # nm      parameter of model
+b = 1.0  # eV      regular
+c = 0.01  # nm
+
+# Helpful solution to convert itertools combinations to numpy arrays here:
+## https://stackoverflow.com/questions/33282369/convert-itertools-array-into-numpy-array
+def cp(l):
+    return np.fromiter(
+        itertools.chain(*itertools.combinations(l, 2)), dtype=int
+    ).reshape(-1, 2)
+
+
+class Cluster:
+    def __init__(self, r_na, r_cl):
+        """
+        Inputs the list of Na and Cl positions. Na has charge +1, Cl has -1.
+        The array of ions itself does not change throughout the calculation, and
+        neither do the charges. As such, we can just compute the combinations one time
+        and refer to it throughout the calculation.
+        """
+        self.positions = np.concatenate((r_na, r_cl))
+        self.charges = np.concatenate(
+            [np.ones(r_na.shape[0]), np.full(r_cl.shape[0], -1)]
+        )
+        self.combs = cp(np.arange(self.charges.size))
+        self.chargeprods = (
+            self.charges[self.combs][:, 0] * self.charges[self.combs][:, 1]
+        )
+        self.rij = np.linalg.norm(
+            self.positions[self.combs][:, 0] - self.positions[self.combs][:, 1], axis=1
+        )
+
+    def Vij(self):
+        """Calculate a numpy vector of all of the potentials of the combinations"""
+        self.Vij_ = np.zeros_like(self.rij)
+        pos = self.chargeprods > 0
+        neg = ~pos
+        self.Vij_[pos] = ke2 / self.rij[pos] + b * (c / self.rij[pos]) ** 12
+        self.Vij_[neg] = (
+            -ke2 / self.rij[neg]
+            + alpha * np.exp(-self.rij[neg] / rho)
+            + b * (c / self.rij[neg]) ** 12
+        )
+        return self.Vij_
+
+    def V(self):
+        """Total potential, which is a sum of the Vij vector"""
+        return np.sum(self.Vij())
+
+    def get_vals(self):
+        """Positions interpreted as a flat shape"""
+        return np.reshape(self.positions, -1)
+
+    def set_vals(self, vals):
+        """Inputs flat shape of positions, used by __call__"""
+        self.positions = vals.reshape(self.positions.shape)
+        self.rij = np.linalg.norm(
+            self.positions[self.combs][:, 0] - self.positions[self.combs][:, 1], axis=1
+        )
+
+    def __call__(self, vals):
+        """Function that  scipy.optimize.minimize will call"""
+        self.set_vals(vals)
+        return self.V()
+
+
+a = 0.2
+
+# r_na = np.array([[0, 0, 0], [0, 0, 2*a],[0,a,a],[0,2*a,-0.5*a]])
+# r_cl = np.array([[0, 0, a], [0, a, 2*a],[a,a,0],[a,- a,0]])
+
+while True:
+
+    r_na = np.random.rand(4, 3)*.75
+    r_cl = np.random.rand(4, 3)*.75
+
+    cluster = Cluster(r_na, r_cl)
+    vals_init = cluster.get_vals()
+    # print("initial Na positions:\n", r_na)
+    # print("initial Cl positions:\n", r_cl)
+    # print("initial positions flattened shape:\n", vals_init)
+    # print("initial V  :", cluster.V())
+
+    res = scipy.optimize.minimize(fun=cluster, x0=vals_init, tol=1e-3, method="BFGS")
+    cluster.set_vals(
+        res.x
+    )  # For some reason, "minimize" is not updating the class at the last iteration
+
+    if (
+        res.fun < -26.9999
+        and not np.isclose(res.fun, -27.7298, 0.0005)  # rectangular potential
+        and not np.isclose(res.fun, -28.2358, 0.0005)  # cubic potential
+        and not np.isclose(res.fun, -26.1369, 0.0005)  # just a straight line...
+        and not np.isclose(res.fun, -27.7824, 0.0005)  # octagon
+        and not np.isclose(res.fun, -27.3422, 0.0005)
+    ):
+        break
+print("Final optimized cluster positions")
+print(cluster.positions)
+print("Final potential:", res.fun)
+
+
+
+fig = plt.figure()
+ax = fig.add_subplot(111, projection="3d")
+
+charges = cluster.charges
+x, y, z = cluster.positions[:, 0], cluster.positions[:, 1], cluster.positions[:, 2]
+ax.scatter(x, y, z, c=charges, cmap="coolwarm", s = 200)
+ax.set_xlabel("x")
+ax.set_ylabel("y")
+ax.set_zlabel("z")
+plt.show()