Removed deprecated code
This commit is contained in:
parent
a559f72221
commit
ba893881a3
@ -57,15 +57,15 @@ def occupation_matrix(
|
||||
|
||||
|
||||
def _calc_histogram(
|
||||
numberlist: ArrayLike, trajectory: Coordinates, bins: ArrayLike
|
||||
indices: ArrayLike, trajectory: Coordinates, bins: ArrayLike
|
||||
) -> NDArray:
|
||||
matbin = None
|
||||
for index in range(0, len(numberlist), 1000):
|
||||
for index in range(0, len(indices), 1000):
|
||||
try:
|
||||
indices = numberlist[index : index + 1000]
|
||||
current_indices = indices[index : index + 1000]
|
||||
except IndexError:
|
||||
indices = numberlist[index:]
|
||||
frames = np.concatenate(np.array([trajectory.pbc[i] for i in indices]))
|
||||
current_indices = indices[index:]
|
||||
frames = np.concatenate(np.array([trajectory.pbc[i] for i in current_indices]))
|
||||
hist, _ = np.histogramdd(frames, bins=bins)
|
||||
if matbin is None:
|
||||
matbin = hist
|
||||
@ -202,436 +202,6 @@ def find_energy_maxima(
|
||||
x = np.array(data_d["r"])
|
||||
y = np.array(data_d["energy"])
|
||||
mask = (x >= r_min) * (x <= r_max)
|
||||
print(x[mask], y[mask])
|
||||
p3 = Poly.fit(x[mask], y[mask], deg=2)
|
||||
print(p3)
|
||||
energies.append(np.max(p3(np.linspace(r_min, r_max, 1000))))
|
||||
return pd.DataFrame({"d": distances, "energy": energies})
|
||||
|
||||
|
||||
# DEPRECATED CODE STARTING FROM HERE
|
||||
# LEFT FOR TESTING ONLY
|
||||
VALID_GEOMETRY = {"cylindrical", "slab"}
|
||||
|
||||
|
||||
def get_fel(
|
||||
traj,
|
||||
path,
|
||||
geometry,
|
||||
temperature,
|
||||
edge=0.05,
|
||||
radiusmin=0.05,
|
||||
radiusmax=2.05,
|
||||
z=[-np.inf, np.inf],
|
||||
overwrite=False,
|
||||
):
|
||||
"""
|
||||
The main method of this script. Will calculate the energy difference based on radius.
|
||||
|
||||
This method will calculate the relative energy of different minima based on radius
|
||||
to the pore center.
|
||||
After that it will save those results in a .npy file in the filepath give by the
|
||||
"path" parameter and will also try to load from there.
|
||||
|
||||
|
||||
Parameters:
|
||||
traj : The trajectory of the system to be evaluated
|
||||
path : The save and load location of the files
|
||||
geometry : Either "cylindrical" or "slab". The geometry of your system. Other types
|
||||
currently not supported
|
||||
temperature: The temperature of your system. Needed for the energy difference
|
||||
edge (opt.) : The length of the cubes in which your system will be divided
|
||||
radiusmin (opt.) : The radius where the calculation begins. Will create a bin of
|
||||
+- 0.05 of that number.
|
||||
radiusmax (opt.) : The radius where the calculation ends. Will create a bin of
|
||||
+- 0.05 of that number.
|
||||
z (opt.) : The evaluated slice of the trajectory for the energy landscape.
|
||||
|
||||
Returns:
|
||||
list: A list of the energy difference based on radius
|
||||
|
||||
"""
|
||||
if geometry not in VALID_GEOMETRY:
|
||||
raise ValueError("results: status must be one of %r." % VALID_GEOMETRY)
|
||||
|
||||
if (os.path.exists(f"{path}/radiiData.npy")) & (not overwrite):
|
||||
data = np.load(f"{path}/radiiData.npy")
|
||||
bins = np.load(f"{path}/radiiBins.npy")
|
||||
# Here the different standard geometries are inserted
|
||||
else:
|
||||
if geometry == "cylindrical":
|
||||
bins, data = short_all_radii(
|
||||
traj, path, edge=edge, radiusmin=radiusmin, radiusmax=radiusmax, z=z
|
||||
)
|
||||
elif geometry == "slab":
|
||||
bins, data = short_all_radii_slab(
|
||||
traj, path, edge=edge, radiusmin=radiusmin, radiusmax=radiusmax, z=z
|
||||
)
|
||||
np.save(f"{path}/radiiData", data)
|
||||
np.save(f"{path}/radiiBins", bins)
|
||||
|
||||
energy_differences = np.array(calculate_energy_difference(data, bins, temperature))
|
||||
d = np.linspace(radiusmin, radiusmax, len(energy_differences))
|
||||
return d, energy_differences
|
||||
|
||||
|
||||
def fill_bins(traj, path, edge=0.05):
|
||||
# If available Matrix is directly loaded
|
||||
if os.path.exists(f"{path}/Matrix{edge}.npy"):
|
||||
matbin = np.load(f"{path}/Matrix{edge}.npy")
|
||||
return matbin
|
||||
|
||||
pool = mp.Pool(8)
|
||||
size = math.ceil(len(traj) / 8)
|
||||
|
||||
# Trajectory is split for parallel computing
|
||||
indices = list(chunksplit(np.arange(0, len(traj), 1), size))
|
||||
# indices = list(Chunksplit(np.arange(len(traj)-80, len(traj), 1), size))
|
||||
fill = partial(help_fill, traj=traj)
|
||||
results = pool.map(fill, indices)
|
||||
boxmat = traj[0].box
|
||||
|
||||
a = math.ceil(boxmat[0][0] / 0.05)
|
||||
b = math.ceil(boxmat[1][1] / 0.05)
|
||||
c = math.ceil(boxmat[2][2] / 0.05)
|
||||
matbin = np.zeros((a, b, c))
|
||||
|
||||
pool.close()
|
||||
|
||||
for mat in results:
|
||||
matbin = matbin + mat
|
||||
np.save(file=f"{path}/Matrix{edge}", arr=matbin)
|
||||
return matbin
|
||||
|
||||
|
||||
def help_fill(numberlist, traj, edge=0.05):
|
||||
boxmat = traj[0].box
|
||||
|
||||
a = math.ceil(boxmat[0][0] / edge)
|
||||
b = math.ceil(boxmat[1][1] / edge)
|
||||
c = math.ceil(boxmat[2][2] / edge)
|
||||
matbin = np.zeros((a, b, c))
|
||||
temp = np.array([[]]).reshape(0, 3)
|
||||
|
||||
# Trajectory is split in chunks of 1000 frames to increase efficency while
|
||||
# keeping ram usage low
|
||||
h = 1000
|
||||
while h < len(numberlist):
|
||||
temp = np.array([[]]).reshape(0, 3)
|
||||
x = numberlist[h - 1000]
|
||||
y = numberlist[h]
|
||||
for j in traj.pbc[x:y]:
|
||||
l = np.floor(j / edge).astype("int32")
|
||||
temp = np.concatenate((temp, np.array(l)))
|
||||
# Positions are counted for whole chunk
|
||||
unique, counts = np.unique(temp, return_counts=True, axis=0)
|
||||
m = 0
|
||||
# Count is combined into matrix with position in Matrix corresponding to
|
||||
# position in system
|
||||
for z in unique.astype("int"):
|
||||
a = z[0]
|
||||
b = z[1]
|
||||
c = z[2]
|
||||
matbin[a][b][c] += counts[m]
|
||||
m += 1
|
||||
h += 1000
|
||||
# The last few frames of the system are seperately calculated
|
||||
x = numberlist[h - 1000]
|
||||
y = numberlist[-1]
|
||||
temp = np.array([[]]).reshape(0, 3)
|
||||
for j in traj.pbc[x : y + 1]:
|
||||
l = np.floor(j / edge).astype("int32")
|
||||
temp = np.concatenate((temp, np.array(l)))
|
||||
|
||||
unique, counts = np.unique(temp, return_counts=True, axis=0)
|
||||
m = 0
|
||||
for z in unique.astype("int"):
|
||||
a = z[0]
|
||||
b = z[1]
|
||||
c = z[2]
|
||||
matbin[a][b][c] += counts[m]
|
||||
m += 1
|
||||
return matbin
|
||||
|
||||
|
||||
def calculate_maxima(matbin):
|
||||
maxima = []
|
||||
|
||||
z = [-1, 0, 1]
|
||||
|
||||
max_i = matbin.shape[0]
|
||||
max_j = matbin.shape[1]
|
||||
max_k = matbin.shape[2]
|
||||
|
||||
# For each element in the matrix all surrounding 26 elements are compared to
|
||||
# determine the minimum.
|
||||
# Algorithm can definitely be improved but is so fast that it's not necessary.
|
||||
|
||||
for i in range(max_i):
|
||||
for j in range(max_j):
|
||||
for k in range(max_k):
|
||||
a = matbin[i][j][k]
|
||||
b = True
|
||||
for l in z:
|
||||
for m in z:
|
||||
for n in z:
|
||||
if (l != 0) or (m != 0) or (n != 0):
|
||||
b = b and (
|
||||
a
|
||||
> matbin[(i + l) % max_i][(j + m) % max_j][
|
||||
(k + n) % max_k
|
||||
]
|
||||
)
|
||||
|
||||
if b:
|
||||
maxima.append([i, j, k])
|
||||
return maxima
|
||||
|
||||
|
||||
def list_of_coordinates(matbin):
|
||||
# Matrix elements are translated back into postion vectors
|
||||
max_i = matbin.shape[0]
|
||||
max_j = matbin.shape[1]
|
||||
max_k = matbin.shape[2]
|
||||
|
||||
coord = np.zeros((max_i, max_j, max_k, 3))
|
||||
|
||||
for i in range(max_i):
|
||||
for j in range(max_j):
|
||||
for k in range(max_k):
|
||||
coord[i][j][k] = [i, j, k]
|
||||
|
||||
return coord
|
||||
|
||||
|
||||
def sphere_quotient(
|
||||
matbin,
|
||||
maxima,
|
||||
coordlist,
|
||||
tree,
|
||||
radius,
|
||||
box,
|
||||
edge,
|
||||
unitdist,
|
||||
z=np.array([-np.inf, np.inf]),
|
||||
):
|
||||
# Here pore center is assumed to be system center
|
||||
diff = np.diag(box)[0:2] / 2
|
||||
|
||||
# Distance to the z-axis
|
||||
distance = np.linalg.norm(np.array(maxima)[:, 0:2] * edge - diff, axis=1)
|
||||
|
||||
# selection with given parameters
|
||||
mask = (
|
||||
(distance > radius - 0.05)
|
||||
& (distance < radius + 0.05)
|
||||
& (np.array(maxima)[:, 2] * edge > z[0])
|
||||
& (np.array(maxima)[:, 2] * edge < z[1])
|
||||
)
|
||||
|
||||
maxima_masked = np.array(maxima)[mask]
|
||||
|
||||
# Distances between maxima and other cubes in the system
|
||||
coordlist = coordlist.reshape(-1, 3)
|
||||
numOfNeigbour = np.max(tree.query_ball_point(maxima, unitdist, return_length=True))
|
||||
d, neighbourlist = tree.query(maxima_masked, k=numOfNeigbour, workers=-1)
|
||||
|
||||
i = 0
|
||||
|
||||
average = 0
|
||||
|
||||
y = []
|
||||
num = []
|
||||
for neighbours in neighbourlist:
|
||||
current = maxima_masked[i]
|
||||
|
||||
# energy between minimum and all sourrounding cubes is calculated
|
||||
energy = -np.log(
|
||||
matbin[current[0], current[1], current[2]]
|
||||
/ matbin.flatten()[neighbours[1:]]
|
||||
)
|
||||
|
||||
v, b = np.histogram(
|
||||
d[i, 1:].flatten()[(energy < np.Inf) & (energy > -np.Inf)],
|
||||
bins=np.arange(1, 40.5, 0.5),
|
||||
weights=energy[(energy < np.Inf) & (energy > -np.Inf)],
|
||||
)
|
||||
|
||||
# For averaging purposes number weighted is also the calculated
|
||||
k, l = np.histogram(
|
||||
d[i, 1:].flatten()[(energy < np.Inf) & (energy > -np.Inf)],
|
||||
bins=np.arange(1, 40.5, 0.5),
|
||||
)
|
||||
y.append(v)
|
||||
num.append(k)
|
||||
i += 1
|
||||
|
||||
# energy is averaged over minima
|
||||
quotient = np.sum(y, axis=0) / np.sum(num, axis=0)
|
||||
if len(neighbourlist) == 0:
|
||||
return np.arange(1, 40.5, 0.5) * edge, np.arange(1, 40, 0.5) * 0
|
||||
return b * edge, quotient
|
||||
|
||||
|
||||
def sphere_quotient_slab(
|
||||
matbin, maxima, coordlist, tree, radius, box, edge, unitdist, z=[-np.inf, np.inf]
|
||||
):
|
||||
# Same as SphereQuotient bur for Slabs
|
||||
diff = box[2, 2] / 2
|
||||
distance = abs(np.array(maxima)[:, 2] * edge - diff)
|
||||
mask = (
|
||||
(distance > radius - 0.05)
|
||||
& (distance < radius + 0.05)
|
||||
& (np.array(maxima)[:, 2] > z[0])
|
||||
& (np.array(maxima)[:, 2] < z[1])
|
||||
)
|
||||
|
||||
maxima_masked = np.array(maxima)[mask]
|
||||
|
||||
coordlist = coordlist.reshape(-1, 3)
|
||||
numOfNeigbour = np.max(tree.query_ball_point(maxima, unitdist, return_length=True))
|
||||
d, neighbourlist = tree.query(maxima_masked, k=numOfNeigbour, workers=-1)
|
||||
|
||||
i = 0
|
||||
|
||||
average = 0
|
||||
|
||||
y = []
|
||||
num = []
|
||||
for neighbours in neighbourlist:
|
||||
current = maxima_masked[i]
|
||||
|
||||
energy = -np.log(
|
||||
matbin[current[0], current[1], current[2]]
|
||||
/ matbin.flatten()[neighbours[1:]]
|
||||
)
|
||||
|
||||
# Energy is ordered according to distance to the minimum
|
||||
v, z = np.histogram(
|
||||
d[i, 1:].flatten()[(energy < np.Inf) & (energy > -np.Inf)],
|
||||
bins=unitdist * 2 - 2,
|
||||
weights=energy[(energy < np.Inf) & (energy > -np.Inf)],
|
||||
)
|
||||
|
||||
k, l = np.histogram(
|
||||
d[i, 1:].flatten()[(energy < np.Inf) & (energy > -np.Inf)],
|
||||
bins=unitdist * 2 - 2,
|
||||
)
|
||||
|
||||
y.append(v)
|
||||
|
||||
num.append(k)
|
||||
i += 1
|
||||
|
||||
quotient = np.sum(y, axis=0) / np.sum(num, axis=0)
|
||||
if len(neighbourlist) == 0:
|
||||
return np.arange(1, 40.5, 0.5), np.arange(1, 40, 0.5) * 0
|
||||
return z * edge, quotient
|
||||
|
||||
|
||||
def short_all_radii(
|
||||
traj, path, edge=0.05, radiusmin=0.05, radiusmax=2.05, z=[-np.inf, np.inf]
|
||||
):
|
||||
# Shorthand function cylindrical systems
|
||||
matbin = fill_bins(traj, path, edge)
|
||||
maxima = calculate_maxima(matbin)
|
||||
coordinates = list_of_coordinates(matbin)
|
||||
tree = KDTree(np.reshape(coordinates, (-1, 3)), boxsize=matbin.shape)
|
||||
bins = []
|
||||
data = []
|
||||
for radius in np.arange(radiusmin, radiusmax, 0.1):
|
||||
b, d = sphere_quotient(
|
||||
matbin,
|
||||
maxima,
|
||||
coordinates,
|
||||
tree,
|
||||
radius=radius,
|
||||
box=traj[0].box,
|
||||
edge=edge,
|
||||
unitdist=40,
|
||||
z=z,
|
||||
)
|
||||
bins.append(b)
|
||||
data.append(d)
|
||||
return bins, data
|
||||
|
||||
|
||||
def short_all_radii_slab(
|
||||
traj, path, edge=0.05, radiusmin=0.05, radiusmax=2.05, z=[-np.inf, np.inf]
|
||||
):
|
||||
# Shorthand function for Slab systems
|
||||
matbin = fill_bins(traj, path, edge)
|
||||
maxima = calculate_maxima(matbin)
|
||||
coordinates = list_of_coordinates(matbin)
|
||||
tree = KDTree(np.reshape(coordinates, (-1, 3)), boxsize=matbin.shape)
|
||||
bins = []
|
||||
data = []
|
||||
for radius in np.arange(radiusmin, radiusmax, 0.1):
|
||||
c, d = sphere_quotient_slab(
|
||||
matbin,
|
||||
maxima,
|
||||
coordinates,
|
||||
tree,
|
||||
radius=radius,
|
||||
box=traj[0].box,
|
||||
edge=edge,
|
||||
unitdist=40,
|
||||
z=z,
|
||||
)
|
||||
bins.append(c)
|
||||
data.append(d)
|
||||
return bins, data
|
||||
|
||||
|
||||
def calculate_energy_difference(data, bins, temperature):
|
||||
"""
|
||||
Calculates the energy difference between local energy minimum and the
|
||||
minimum in the center
|
||||
"""
|
||||
difference = []
|
||||
i = 0
|
||||
while i < len(data):
|
||||
#
|
||||
q = (
|
||||
get_minimum(data[0], bins[0])[1] - get_minimum(data[i], bins[0])[1]
|
||||
) * temperature
|
||||
difference.append(q)
|
||||
i += 1
|
||||
return difference
|
||||
|
||||
|
||||
def chunksplit(list_a, chunk_size):
|
||||
for i in range(0, len(list_a), chunk_size):
|
||||
yield list_a[i : i + chunk_size]
|
||||
|
||||
|
||||
def get_minimum(data, bins):
|
||||
# Fits a polynom of order 3 to determine the energy minimum and analytically
|
||||
# calculates the minimum location
|
||||
y = data[:10]
|
||||
x = bins[1:11]
|
||||
|
||||
popt, _ = scipy.optimize.curve_fit(pol3, x, y, maxfev=80000)
|
||||
|
||||
a, b = solve_quadratic(3 * popt[0], 2 * popt[1], popt[2])
|
||||
|
||||
if 6 * popt[0] * a + 2 * popt[1] > 0:
|
||||
if (np.real(a) < 0) or (np.real(a) > 0.3):
|
||||
return np.real(a), np.average(y)
|
||||
return np.real(a), np.real(pol3(a, *popt))
|
||||
else:
|
||||
if (np.real(b) < 0) or (np.real(b) > 0.3):
|
||||
return np.real(b), np.average(y)
|
||||
return np.real(b), np.real(pol3(b, *popt))
|
||||
|
||||
|
||||
def solve_quadratic(a, b, c):
|
||||
d = (b**2) - (4 * a * c)
|
||||
|
||||
sol1 = (-b - cmath.sqrt(d)) / (2 * a)
|
||||
sol2 = (-b + cmath.sqrt(d)) / (2 * a)
|
||||
return sol1, sol2
|
||||
|
||||
|
||||
def pol3(x, a, b, c, d):
|
||||
return a * x**3 + b * x**2 + c * x + d
|
||||
|
||||
Loading…
Reference in New Issue
Block a user