diff --git a/pymuonsuite/io/output.py b/pymuonsuite/io/output.py index c855ac9..bbf97a3 100644 --- a/pymuonsuite/io/output.py +++ b/pymuonsuite/io/output.py @@ -16,27 +16,6 @@ from soprano.utils import silence_stdio -def write_tensors(tensors, filename, symbols): - """ - Write out a set of 2 dimensional tensors for every atom in a system. - - | Args: - | tensors(Numpy float array, shape: (Atoms, :, :): A list of tensors - | for each atom. - | filename(str): Filename for file. - | symbols(str array): List containing chemical symbol of each atom in - | system. - | - | Returns: Nothing - """ - tensfile = open(filename, "w") - for i in range(np.size(tensors, 0)): - tensfile.write("{0} {1}\n".format(symbols[i], i)) - tensfile.write( - "\n".join(["\t".join([str(x) for x in ln]) for ln in tensors[i]]) + "\n" - ) - - def write_cluster_report(args, params, clusters): if params["clustering_method"] == "hier": diff --git a/pymuonsuite/quantum/vibrational/grid.py b/pymuonsuite/quantum/vibrational/grid.py deleted file mode 100644 index 4723646..0000000 --- a/pymuonsuite/quantum/vibrational/grid.py +++ /dev/null @@ -1,177 +0,0 @@ -""" -Author: Simone Sturniolo and Adam Laverack -""" - - -import random - -import numpy as np -from soprano.collection.generate import linspaceGen - - -# def calc_wavefunction(sigmas, grid_n, sigmaN=3): -# """ -# Calculate harmonic oscillator wavefunction - -# | Args: -# | sigmas(Numpy float array, shape:(axes)): Displacement amplitude along -# | each axis -# | grid_n(int): Number of grid points along each axis -# | sigmaN(int): Number of standard deviations to cover in range -# | (default is 3) -# | Returns: -# | prob_dens (Numpy float array, shape:(grid_n*3)): Probability density of -# | harmonic oscillator at each displacement -# """ -# R_axes = np.array([np.linspace(-sigmas * Ri, sigmas * Ri, grid_n) for Ri in R]) - -# # Wavefunction -# psi_norm = (1.0 / (np.prod(R) ** 2 * np.pi ** 3)) ** 0.25 -# # And along the three axes -# psi = psi_norm * np.exp(-((R_axes / R[:, None]) ** 2) / 2.0) -# # And average -# r2psi2 = R_axes ** 2 * np.abs(psi) ** 2 - -# # Convert to portable output format -# prob_dens = np.zeros((grid_n * 3)) -# for i, mode in enumerate(r2psi2): -# for j, point in enumerate(mode): -# prob_dens[j + i * grid_n] = point - -# return prob_dens - - -def create_displaced_cell(cell, displacements): - """ - Take a set of atomic displacements and a cell and return an ASE Atoms object - with atoms displaced appropriately. - - | Args: - | cell(ASE Atoms object): Cell containing original atomic positions. - | displacements(Numpy float array(num_atoms, 3)): Array containing a - | displacement vector for each atom in the system. - | - | Returns: - | disp_cell(ASE Atoms object): Cell containing displaced atomic positions. - """ - disp_pos = cell.get_positions() - disp_cell = cell.copy() - disp_pos += displacements - disp_cell.set_positions(disp_pos) - disp_cell.calc = cell.calc - - return disp_cell - - -def displaced_cell_range(cell, a_i, grid_n, disp): - """Return a generator of ASE Atoms objects with the displacement of the atom - at index a_i varying between -disp and +disp with grid_n increments - - | Args: - | cell (ASE Atoms object): Object containing atom to be displaced - | a_i (int): Index of atom to be displaced - | grid_n (int): Number of increments/objects to create - | disp (float): Maximum displacement from original position - | - | Returns: - | lg(Soprano linspaceGen object): Generator of displaced cells - """ - pos = cell.get_positions() - cell_L = cell.copy() - pos_L = pos.copy() - pos_L[a_i] -= disp - cell_L.set_positions(pos_L) - cell_R = cell.copy() - pos_R = pos.copy() - pos_R[a_i] += disp - cell_R.set_positions(pos_R) - lg = linspaceGen(cell_L, cell_R, steps=grid_n, periodic=True) - return lg - - -def tl_disp_generator(norm_coords, evecs, num_atoms): - """ - Calculate a set of displacements of atoms in a system by generating a set of - random thermal lines at T=0. - - | Args: - | norm_coords(Numpy float array(num_atoms, num_modes)): Array containing - | the normal mode coordinates of each real mode of the system. - | evecs(Numpy float array(size(norm_coords), num_atoms)): Array containing - | the eigenvectors of all real phonon modes for each atom in the - | system in the format evecs[modes][atoms]. - | num_atoms(int): Number of atoms in the system. - | - | Returns: - | displacements(Numpy float array(num_atoms, 3)): Array containing the - | appropriate displacements of atoms for a randomly generated thermal - | line. - """ - displacements = np.zeros((num_atoms, 3)) - coefficients = np.zeros(np.size(norm_coords, 1)) - for i in range(np.size(coefficients)): - coefficients[i] = random.choice([-1, 1]) - norm_coords = norm_coords * coefficients - - for atom in range(num_atoms): - for mode in range(np.size(norm_coords, 1)): - displacements[atom] += ( - norm_coords[atom][mode] * evecs[mode][atom].real * 1e10 - ) - - return displacements - - -def weighted_tens_avg(tensors, weight): - """ - Given a set of tensors resulting from the sampling of a property over a - set of different displacements, calculate a weighted average of the tensors - for each atom or the whole system, using a given weight for each - displacement. - - | Args: - | tensors(Numpy float array, shape:(N,No. of atoms,x,y)): For each grid - | point, a set of x by y tensors for each atom or the whole system. - | "No. of atoms" should be 1 in the latter case. - | weight(Numpy float array, shape:(N)): A weighting for each point - | on the grid. - | - | Returns: - | tens_avg(Numpy float array, shape:(Atoms,x,y)): The averaged tensor for - | each atom. - """ - num_atoms = np.size(tensors, 1) - x = np.size(tensors, 2) - y = np.size(tensors, 3) - tens_avg = np.zeros((num_atoms, x, y)) - tensors = tensors * weight[:, None, None, None] - for i in range(num_atoms): - tens_avg[i] = np.sum(tensors[:, i], axis=0) / np.sum(weight) - return tens_avg - - -def wf_disp_generator(disp_factor, maj_evecs, grid_n): - """ - Generate a set of displacements of an atom for the wavefunction sampling - method. - - | Args: - | disp_factor(float(3)): A displacement factor for each of the 3 major - | phonon modes of the atom. - | maj_evecs(Numpy float array(3, 3)): The eigenvectors of the 3 major - | phonon modes of the atom. - | grid_n(int): The number of desired grid points along each mode. - | - | Returns: - | displacements(Numpy float array(grid_n*3, 3)): Array containing the - | amount the atom should be displaced by at each grid point. The first - | elements are for the first mode, the second for - | the second mode, etc. - """ - displacements = np.zeros((grid_n * 3, 3)) - max_disp = 3 * maj_evecs * disp_factor[:, None] - for mode in range(3): - for n, t in enumerate(np.linspace(-1, 1, grid_n)): - displacements[n + mode * grid_n] = t * max_disp[mode] - - return displacements diff --git a/pymuonsuite/quantum/vibrational/reports.py b/pymuonsuite/quantum/vibrational/reports.py deleted file mode 100644 index 7c34100..0000000 --- a/pymuonsuite/quantum/vibrational/reports.py +++ /dev/null @@ -1,90 +0,0 @@ -""" -Author: Adam Laverack and Simone Sturniolo -""" - - -import numpy as np -import scipy.constants as cnst - - -def harm_potential_report(R, grid_n, mass, freqs, E_table, filename): - """ - Calculate the harmonic potential at all displacements on the grid for an - atom and write out to file in a format that can be plotted. - - | Args: - | R(Numpy float array, shape:(axes)): Displacement amplitude along each - | axis - | grid_n(int): Number of grid points along each axis - | mass(float): Mass of atom - | freqs(Numpy float array, shape:(axes)): Frequencies of harmonic - | oscillator along each axis - | E_table(Numpy float array, shape:(np.size(R), grid_n)): Table of CASTEP - | final system energies. - | filename(str): Filename to be used for file - | - | Returns: Nothing - """ - R_axes = np.array([np.linspace(-3 * Ri, 3 * Ri, grid_n) for Ri in R]) - # Now the potential, measured vs. theoretical - harm_K = mass * freqs**2 - harm_V = (0.5 * harm_K[:, None] * (R_axes * 1e-10) ** 2) / cnst.electron_volt - # Normalise E_table - if E_table.shape[1] % 2 == 1: - E_table -= (E_table[:, E_table.shape[1] // 2])[:, None] - else: - E_table -= ( - E_table[:, E_table.shape[1] // 2] + E_table[:, E_table.shape[1] // 2 - 1] - )[:, None] / 2.0 - all_table = np.concatenate((R_axes, harm_V, E_table), axis=0) - np.savetxt(filename, all_table.T) - - -def hfine_report(total_grid_n, tensors, hfine_tens_avg, weight, filename, atoms): - """Write a report on a selection of atom's hyperfine coupling constants and - their hyperfine tensor dipolar components based on a vibrational averaging - calculation. - - | Args: - | total_grid_n(int): Total number of grid points - | tensors(Numpy float array, shape:(total_grid_n, num_atoms, 3, 3)): - | Array of hyperfine tensors for each atom at each grid point - | hfine_tens_avg(Numpy float array, shape:(num_atoms,3,3)): Average - | tensors of atoms over grid - | weight (Numpy float, shape:(total_grid_n)): Weighting of each grid - | point - | filename(str): Filename to be used for file - | atoms(dict, {index(int):symbol(str)}): Dictionary containing indices - | and symbols of atoms to write - | hyperfine coupling report about - | - | Returns: - | Nothing - """ - ofile = open(filename, "w") - for index in atoms: - hfine_table = np.trace(tensors[:, index], axis1=1, axis2=2) / 3 - hfine_avg = np.sum(weight * hfine_table) / np.sum(weight) - ofile.write( - "Predicted hyperfine coupling on labeled atom ({1}): {0} MHz\n".format( - hfine_avg, atoms[index] - ) - ) - - evals, evecs = np.linalg.eigh(hfine_tens_avg[index]) - evals, evecs = zip(*sorted(zip(evals, evecs), key=lambda x: abs(x[0]))) - evals_notr = -np.array(evals) + np.average(evals) - - if abs(evals_notr[2]) > abs(evals_notr[0]): - D1 = evals_notr[2] - D2 = evals_notr[1] - evals_notr[0] - else: - D1 = evals_notr[0] - D2 = evals_notr[2] - evals_notr[1] - - ofile.write( - ( - "Predicted dipolar hyperfine components on labeled atom ({2}):\n" - "D1:\t{0} MHz\nD2:\t{1} MHz\n" - ).format(D1, D2, atoms[index]) - ) diff --git a/pymuonsuite/utils.py b/pymuonsuite/utils.py index a5513b3..470333d 100644 --- a/pymuonsuite/utils.py +++ b/pymuonsuite/utils.py @@ -79,48 +79,6 @@ def make_process_slices(N, M): return slices -def create_plane_grid(hkl, cell, f0, f1, N=20): - - # Create a grid of points along a crystal plane - - hkl = np.array(hkl).astype(int) - f0 = np.array(f0) - f1 = np.array(f1) - - # First: verify that the given points (given in fractional coordinates) - # DO belong to the same plane - f01 = f1 - f0 - if np.isclose(np.linalg.norm(f01), 0): - raise ValueError("Points f0 and f1 are too close") - if not np.isclose(np.dot(hkl, f01), 0): - raise ValueError("Points f0 and f1 do not belong to the same plane") - - # Now move to direct space - n = np.dot(hkl, np.linalg.inv(cell)) - p0 = np.dot(cell.T, f0) - p1 = np.dot(cell.T, f1) - - n /= np.linalg.norm(n) - - # Find the scanning directions - p01 = p1 - p0 - - plx = np.zeros(3) - plx[np.where(p01 != 0)[0][0]] = 1 - ply = p01 - np.dot(p01, plx) * plx - ply /= np.linalg.norm(ply) - ply *= np.sign(ply[np.where(ply != 0)[0][0]]) - - # Now to actually create the scanning grid - plgrid = np.array(np.meshgrid(*[np.linspace(0, 1, N)] * 2, indexing="ij")) - - xyzgrid = plgrid[0, None] * plx[:, None, None] * np.dot(p01, plx) - xyzgrid += plgrid[1, None] * ply[:, None, None] * np.dot(p01, ply) - xyzgrid += p0[:, None, None] - - return xyzgrid - - class BackupFile: """Backup a file before performing an operation