matscipy package¶

Submodules¶

matscipy.angle_distribution module¶

matscipy.angle_distribution.angle_distribution(i, j, dr, nbins, *args)¶

Compute a bond angle distribution from a neighbour list.

Parameters

j, dr (i,) – Neighbour list, including list of distance vectors.
nbins (int) – Number of bins for bond angle histogram.
cutoff (float, optional) – Bond length cutoff, i.e. consider only bonds shorter than this length.

matscipy.atomic_strain module¶

Compute deformation gradient tensor and D^2_min measure for non-affine displacements. See: Falk, Langer, Phys. Rev. E 57, 7192 (1998)

matscipy.atomic_strain.array_inverse(A)¶: Compute inverse for each matrix in a list of matrices. This is faster than calling numpy.linalg.inv for each matrix.

matscipy.atomic_strain.atomic_strain(atoms_now, atoms_old, cutoff=None, neighbours=None)¶

Calculate deformation gradient tensor and D^2_min measure for non-affine displacements. See: Falk, Langer, Phys. Rev. B 57, 7192 (1998)

atoms_nowase.Atoms: Current atomic configuration
atoms_oldase.Atoms: Reference atomic configuration
cutofffloat: Neighbor list cutoff.
neighbours( array_like, array_like ): Neighbor list. Automatically computed if not provided.

delta_plus_epsilonarray: 3x3 deformation gradient tensor for each atom.
residualarray: D^2_min norm for each atom

matscipy.atomic_strain.get_D_square_min(atoms_now, atoms_old, i_now, j_now, delta_plus_epsilon=None)¶: Calculate the D^2_min norm of Falk and Langer

matscipy.atomic_strain.get_XIJ(nat, i_now, dr_now, dr_old)¶: Calculates the X_{ij} matrix

matscipy.atomic_strain.get_YIJ(nat, i_now, dr_old)¶: Calculates the Y_{ij} matrix

matscipy.atomic_strain.get_delta_plus_epsilon(nat, i_now, dr_now, dr_old)¶: Calculate delta_ij+epsilon_ij, i.e. the deformation gradient matrix

matscipy.atomic_strain.get_delta_plus_epsilon_dgesv(nat, i_now, dr_now, dr_old)¶: Calculate delta_ij+epsilon_ij, i.e. the deformation gradient matrix

matscipy.elasticity module¶

class matscipy.elasticity.CubicElasticModuli(C11, C12, C44)¶

Bases: object

compliance()¶: Return the compliance coefficients

rotate(A)¶: Compute the rotated stiffness matrix

stiffness()¶: Return the elastic constants

tol = 1e-06¶

matscipy.elasticity.Voigt_6_to_full_3x3_strain(strain_vector)¶: Form a 3x3 strain matrix from a 6 component vector in Voigt notation

matscipy.elasticity.Voigt_6_to_full_3x3_stress(stress_vector)¶: Form a 3x3 stress matrix from a 6 component vector in Voigt notation

matscipy.elasticity.Voigt_6x6_to_cubic(C)¶: Convert the Voigt 6x6 representation into the cubic elastic constants C11, C12 and C44.

matscipy.elasticity.Voigt_6x6_to_full_3x3x3x3(C)¶

Convert from the Voigt representation of the stiffness matrix to the full 3x3x3x3 representation.

Parameters: C (array_like) – 6x6 stiffness matrix (Voigt notation).
Returns: C – 3x3x3x3 stiffness matrix.
Return type: array_like

matscipy.elasticity.cubic_to_Voigt_6x6(C11, C12, C44)¶

matscipy.elasticity.elastic_moduli(C, l=array([1, 0, 0]), R=None, tol=1e-06)¶

Calculate elastic moduli from 6x6 elastic constant matrix C_{ij}.

The elastic moduli calculated are: Young’s muduli, Poisson’s ratios, shear moduli, bulk mudulus and bulk mudulus tensor.

If a direction l is specified, the system is rotated to have it as its x direction (see Notes for details). If R is specified the system is rotated according to it.

Parameters

C (array_like) – 6x6 stiffness matrix (Voigt notation).
l (array_like, optional) – 3D direction vector for pull (the default is the x direction of the original system)
R (array_like, optional) – 3x3 rotation matrix.

Returns

E (array_like) – Young’s modulus for a stress in each of the three directions of the rotated system.
nu (array_like) – 3x3 matrix with Poisson’s ratios.
Gm (array_like) – 3x3 matrix with shear moduli.
B (float) – Bulk modulus.
K (array_like) – 3x3 matrix with bulk modulus tensor.

Other Parameters

tol (float, optional) – tolerance for checking validity of rotation and comparison of vectors.

Notes

It works by rotating the elastic constant tensor to the desired direction, so it should be valid for arbitrary crystal structures. If only l is specified there is an infinite number of possible rotations. The chosen one is a rotation along the axis orthogonal to the plane defined by the vectors (1, 0, 0) and l.

Bulk modulus tensor as defined in O. Rand and V. Rovenski, “Analytical Methods in Anisotropic Elasticity”, Birkh”auser (2005), pp. 71.

matscipy.elasticity.fit_elastic_constants(a, symmetry='triclinic', N_steps=5, delta=0.01, optimizer=None, verbose=True, graphics=False, logfile=None, **kwargs)¶

Compute elastic constants by linear regression of stress vs. strain

Parameters

a (ase.Atoms or list of ase.Atoms) – Either a single atomic configuration or a list of strained configurations. If a single configuration is given, it is passed generate_strained_configs() along with symmetry, N_steps, and delta to generate the set of strained configurations.
symmetry (string) – Symmetry to use to determine which strain patterns are necessary. Default is ‘triclininc’, i.e. no symmetry.
N_steps (int) – Number of atomic configurations to generate for each strain pattern. Default is 5. Absolute strain values range from -delta*N_steps/2 to +delta*N_steps/2.
delta (float) – Strain increment for analytical derivatives of stresses. Default is 1e-2.
optimizer (ase.optimizer.*) – Optimizer to use for atomic positions (internal degrees of freedom) for each applied strain. Initial config a is not optimised, and should already have relaxed cell and atomic positions. Does not optimize atomic positions if set to None.
verbose (bool) – If True, print additional infomation about the quality of fit and summarise results of C_ij and estimated errors. Default True.
graphics (bool) – If True, use matplotlib.pyplot to plot the stress vs. strain curve for each C_ij component fitted. Default True.
logfile (bool) – Log file to write optimizer output to. Default None (i.e. suppress output).
**kwargs (dict) – Additional arguments to pass to optimizer.run() method e.g. fmax.

Returns

C (array_like) – 6x6 matrix of the elastic constants in Voigt notation.
C_err (array_like) – If scipy.stats module is available then error estimates for each C_ij component are obtained from the accuracy of the linear regression. Otherwise an array of np.zeros((6,6)) is returned.

Notes

Code originally adapted from elastics.py script, available from http://github.com/djw/elastic-constants

matscipy.elasticity.full_3x3_to_Voigt_6_index(i, j)¶

matscipy.elasticity.full_3x3_to_Voigt_6_strain(strain_matrix)¶: Form a 6 component strain vector in Voigt notation from a 3x3 matrix

matscipy.elasticity.full_3x3_to_Voigt_6_stress(stress_matrix)¶: Form a 6 component stress vector in Voigt notation from a 3x3 matrix

matscipy.elasticity.full_3x3x3x3_to_Voigt_6x6(C)¶: Convert from the full 3x3x3x3 representation of the stiffness matrix to the representation in Voigt notation. Checks symmetry in that process.

matscipy.elasticity.generate_strained_configs(at0, symmetry='triclinic', N_steps=5, delta=0.01)¶

Generate a sequence of strained configurations

Parameters

a (ase.Atoms) – Bulk crystal configuration - both unit cell and atomic positions should be relaxed before calling this routine.
symmetry (string) – Symmetry to use to determine which strain patterns are necessary. Default is ‘triclininc’, i.e. no symmetry.
N_steps (int) – Number of atomic configurations to generate for each strain pattern. Default is 5. Absolute strain values range from -delta*N_steps/2 to +delta*N_steps/2.
delta (float) – Strain increment for analytical derivatives of stresses. Default 1e-2

Returns

Generator which yields a sequence of ase.Atoms instances corresponding
to the minima set of strained conifugurations required to evaluate the
full 6x6 C_ij matrix under the assumed symmetry.

matscipy.elasticity.invariants(s, syy=None, szz=None, syz=None, sxz=None, sxy=None, full_3x3_to_Voigt_6=<function full_3x3_to_Voigt_6_stress>)¶

matscipy.elasticity.measure_triclinic_elastic_constants(a, delta=0.001, optimizer=None, logfile=None, **kwargs)¶

Brute-force measurement of elastic constants for a triclinic (general) unit cell.

Parameters

a (ase.Atoms) – Atomic configuration.
optimizer (ase.optimizer.*) – Optimizer to use for atomic position. Does not optimize atomic position if set to None.
delta (float) – Strain increment for analytical derivatives of stresses.

Returns

C – 6x6 matrix of the elastic constants in Voigt notation.

Return type

array_like

matscipy.elasticity.poisson_ratio(C, l, m)¶: Calculate approximate Poisson ratio

u_{lm} from 6x6 elastic constant matrix C_{ij}

This is the response in m direction to pulling in l direction. Result is dimensionless.

Formula is from W. Brantley, Calculated elastic constants for stress problems associated with semiconductor devices. J. Appl. Phys., 44, 534 (1973).

matscipy.elasticity.rotate_cubic_elastic_constants(C11, C12, C44, A, tol=1e-06)¶

Return rotated elastic moduli for a cubic crystal given the elastic constant in standard C11, C12, C44 notation.

Parameters

C12, C44 (C11,) – Cubic elastic moduli.
A (array_like) – 3x3 rotation matrix.

Returns

C – 6x6 matrix of rotated elastic constants (Voigt notation).

Return type

array

matscipy.elasticity.rotate_elastic_constants(C, A, tol=1e-06)¶

Return rotated elastic moduli for a general crystal given the elastic constant in Voigt notation.

Parameters

C (array_like) – 6x6 matrix of elastic constants (Voigt notation).
A (array_like) – 3x3 rotation matrix.

Returns

C – 6x6 matrix of rotated elastic constants (Voigt notation).

Return type

array

matscipy.elasticity.youngs_modulus(C, l)¶

Calculate approximate Youngs modulus E_l from 6x6 elastic constants matrix C_ij

This is the modulus for loading in the l direction. For the exact answer, taking into account elastic anisotropuy, rotate the C_ij matrix to the correct frame, compute the compliance matrix:

C = ...  # 6x6 C_ij matrix in crystal frame
A = ...  # rotation matrix
Cr = rotate_elastic_constants(C, A)
S = np.inv(Cr)
E_x = 1/S[0, 0]  # Young's modulus for a pull in x direction
E_y = 1/S[1, 1]  # Young's modulus for a pull in y direction
E_z = 1/S[0, 0]  # Young's modulus for a pull in z direction

Notes

Formula is from W. Brantley, Calculated elastic constants for stress problems associated with semiconductor devices. J. Appl. Phys., 44, 534 (1973).

matscipy.hydrogenate module¶

matscipy.hydrogenate.hydrogenate(a, cutoff, bond_length, b=None, mask=[True, True, True], exclude=None, vacuum=None)¶

Hydrogenate a slab of material at its periodic boundary conditions. Boundary conditions are turned into nonperiodic.

Parameters

a (ase.Atoms) – Atomic configuration.
cutoff (float) – Cutoff for neighbor counting.
bond_length (float) – X-H bond length for hydrogenation.
b (ase.Atoms, optional) – If present, this is the configuration to hydrogenate. Number of atoms must be identical to a object. All bonds present in a but not present in b will be hydrogenated in b.
mask (list of bool) – Cartesian directions which to hydrogenate, only if b argument is not given.
exclude (array_like) – Boolean array masking atoms to be excluded from hydrogenation.
vacuum (float, optional) – Add this much vacuum after hydrogenation.

Returns

a – Atomic configuration of the hydrogenated slab.

Return type

ase.Atoms

matscipy.io module¶

matscipy.logger module¶

Log status to screen.

class matscipy.logger.Logger(logfile=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='UTF-8'>, outevery=1, sepevery=10)¶

Bases: object

flush()¶

get_logfile()¶

has_logfile()¶

iteration_finished()¶

pr(s, caller=None, logfile=None)¶

set_logfile(logfile)¶

set_outevery(outevery)¶

st(hdr, vals, force_print=False)¶

warn(s, caller=None)¶

matscipy.logger.flatten(x)¶

matscipy.logger.hdr_str(s, x)¶: Return header description strings

matscipy.logger.hdrfmt_str(x, i)¶: Return header format string for datatype x

matscipy.logger.numfmt_str(x, i)¶: Return numeric format string for datatype x

matscipy.neighbours module¶

matscipy.neighbours.find_indices_of_reversed_pairs(i_n, j_n, abs_dr_n)¶

Find neighbor list indices where reversed pairs are stored

Given list of identifiers of neighbor atoms i_n and j_n, determines the list of indices reverse into the neighbor list where each pair is reversed, i.e. i_n[reverse[n]]=j_n[n] and j_n[reverse[n]]=i_n[n] for each index n in the neighbor list

In the case of small periodic systems, one needs to be careful, because the same pair may appear more than one time, with different pair distances. Therefore, the pair distance must be taken into account.

We assume that there is in fact one reversed pair for every pair. However, we do not check this assumption in order to avoid overhead.

Parameters

i_n (array_like) – array of atom identifiers
j_n (array_like) – array of atom identifiers
abs_dr_n (array_like) – pair distances

Returns

reverse – array of indices into i_n and j_n

Return type

numpy.ndarray

matscipy.neighbours.mic(dr, cell, pbc=None)¶

Apply minimum image convention to an array of distance vectors.

Parameters

dr (array_like) – Array of distance vectors.
cell (array_like) – Simulation cell.
pbc (array_like, optional) – Periodic boundary conditions in x-, y- and z-direction. Default is to assume periodic boundaries in all directions.

Returns

dr – Array of distance vectors, wrapped according to the minimum image convention.

Return type

array

matscipy.neighbours.neighbour_list(quantities, atoms=None, cutoff=None, positions=None, cell=None, pbc=None, numbers=None, cell_origin=None)¶

Compute a neighbor list for an atomic configuration. Atoms outside periodic boundaries are mapped into the box. Atoms outside nonperiodic boundaries are included in the neighbor list but the complexity of neighbor list search for those can become n^2.

The neighbor list is sorted by first atom index ‘i’, but not by second atom index ‘j’.

The neighbour list accepts either an ASE Atoms object or positions and cell vectors individually.

quantitiesstr
Quantities to compute by the neighbor list algorithm. Each character in this string defines a quantity. They are returned in a tuple of the same order. Possible quantities are

‘i’ : first atom index ‘j’ : second atom index ‘d’ : absolute distance ‘D’ : distance vector ‘S’ : shift vector (number of cell boundaries crossed by the bond

between atom i and j). With the shift vector S, the distances D between atoms can be computed from: D = a.positions[j]-a.positions[i]+S.dot(a.cell)

atomsase.Atoms
Atomic configuration. (Default: None)

cutofffloat or dict

Cutoff for neighbor search. It can be

A single float: This is a global cutoff for all elements.

A dictionary: This specifies cutoff values for element pairs. Specification accepts element numbers of symbols. Example: {(1, 6): 1.1, (1, 1): 1.0, (‘C’, ‘C’): 1.85}

A list/array with a per atom value: This specifies the radius of an atomic sphere for each atoms. If spheres overlap, atoms are within each others neighborhood.

positionsarray_like
Atomic positions. (Default: None)

cellarray_like
Cell vectors as a 3x3 matrix. (Default: Shrink wrapped cell)

pbcarray_like
3-vector containing periodic boundary conditions in all three directions. (Default: Nonperiodic box)

numbersarray_like
Array containing the atomic numbers.

i, j, …array
Tuple with arrays for each quantity specified above. Indices in i are returned in ascending order 0..len(a), but the order of (i,j) pairs is not guaranteed.

Examples assume Atoms object a and numpy imported as np. 1. Coordination counting:

i = neighbor_list(‘i’, a, 1.85) coord = np.bincount(i)

Coordination counting with different cutoffs for each pair of species

i = neighbor_list(‘i’, a,
{(‘H’, ‘H’): 1.1, (‘C’, ‘H’): 1.3, (‘C’, ‘C’): 1.85})

coord = np.bincount(i)

Pair distribution function:
d = neighbor_list(‘d’, a, 10.00) h, bin_edges = np.histogram(d, bins=100) pdf = h/(4*np.pi/3*(bin_edges[1:]**3 - bin_edges[:-1]**3)) * a.get_volume()/len(a)

Pair potential:
i, j, d, D = neighbor_list(‘ijdD’, a, 5.0) energy = (-C/d**6).sum() pair_forces = (6*C/d**5 * (D/d).T).T forces_x = np.bincount(j, weights=pair_forces[:, 0], minlength=len(a)) - np.bincount(i, weights=pair_forces[:, 0], minlength=len(a)) forces_y = np.bincount(j, weights=pair_forces[:, 1], minlength=len(a)) - np.bincount(i, weights=pair_forces[:, 1], minlength=len(a)) forces_z = np.bincount(j, weights=pair_forces[:, 2], minlength=len(a)) - np.bincount(i, weights=pair_forces[:, 2], minlength=len(a))

Dynamical matrix for a pair potential stored in a block sparse format:
from scipy.sparse import bsr_matrix i, j, dr, abs_dr = neighbor_list(‘ijDd’, atoms) energy = (dr.T / abs_dr).T dynmat = -(dde * (energy.reshape(-1, 3, 1) * energy.reshape(-1, 1, 3)).T).T -(de / abs_dr * (np.eye(3, dtype=energy.dtype) - (energy.reshape(-1, 3, 1) * energy.reshape(-1, 1, 3))).T).T dynmat_bsr = bsr_matrix((dynmat, j, first_i), shape=(3*len(a), 3*len(a)))

dynmat_diag = np.empty((len(a), 3, 3)) for x in range(3):

for y in range(3):
dynmat_diag[:, x, y] = -np.bincount(i, weights=dynmat[:, x, y])

dynmat_bsr += bsr_matrix((dynmat_diag, np.arange(len(a)),

np.arange(len(a) + 1)),

shape=(3 * len(a), 3 * len(a)))

i_n, j_n, dr_nc, abs_dr_n = neighbour_list(‘ijDd’, atoms, dict)

e_nc = (dr_nc.T/abs_dr_n).T

D_ncc = -(dde_n * (e_nc.reshape(-1,3,1) * e_nc.reshape(-1,1,3)).T).T D_ncc += -(de_n/abs_dr_n * (np.eye(3, dtype=e_nc.dtype) - (e_nc.reshape(-1,3,1) * e_nc.reshape(-1,1,3))).T).T

D = bsr_matrix((D_ncc, j_n, first_i), shape=(3*nat,3*nat))

Ddiag_icc = np.empty((nat,3,3)) for x in range(3):

for y in range(3):
Ddiag_icc[:,x,y] = -np.bincount(i_n, weights = D_ncc[:,x,y])

D += bsr_matrix((Ddiag_icc,np.arange(nat),np.arange(nat+1)), shape=(3*nat,3*nat))

return D

matscipy.rings module¶

matscipy.rings.ring_statistics(a, cutoff, maxlength=- 1)¶

Compute number of shortest path rings in sample. See: D.S. Franzblau, Phys. Rev. B 44, 4925 (1991)

Parameters

a (ase.Atoms) – Atomic configuration.
cutoff (float) – Cutoff for neighbor counting.
maxlength (float, optional) – Maximum ring length. Search for rings will stop at this length. This is useful to speed up calculations for large systems.

Returns

ringstat – Array with number of shortest path rings.

Return type

array

matscipy.socketcalc module¶

class matscipy.socketcalc.AtomsRequestHandler(request, client_address, server)¶

Bases: socketserver.StreamRequestHandler

handle()¶

matscipy.socketcalc.AtomsServer¶: alias of matscipy.socketcalc.AtomsServerAsync

class matscipy.socketcalc.AtomsServerAsync(server_address, RequestHandlerClass, clients, bind_and_activate=True, max_attempts=3, bgq=False, logger=<matscipy.logger.Logger object>)¶

Bases: matscipy.socketcalc.AtomsServerSync, socketserver.ThreadingMixIn

Asynchronous (threaded) version of AtomsServer

shutdown()¶

Stops the serve_forever loop.

Blocks until the loop has finished. This must be called while serve_forever() is running in another thread, or it will deadlock.

shutdown_clients()¶

class matscipy.socketcalc.AtomsServerSync(server_address, RequestHandlerClass, clients, bind_and_activate=True, max_attempts=3, bgq=False, logger=<matscipy.logger.Logger object>)¶

Bases: socketserver.TCPServer

allow_reuse_address = True¶

get_results()¶

join_all()¶

put(at, client_id, label, force_restart=False)¶

server_activate()¶

Called by constructor to activate the server.

May be overridden.

shutdown()¶

Stops the serve_forever loop.

Blocks until the loop has finished. This must be called while serve_forever() is running in another thread, or it will deadlock.

shutdown_clients()¶

class matscipy.socketcalc.CastepClient(client_id, exe, env=None, npj=1, ppn=1, block=None, corner=None, shape=None, jobname='socketcalc', rundir=None, fmt='REFTRAJ', parmode=None, mpirun='mpirun', mpirun_args=['-np'], logger=<matscipy.logger.Logger object>, max_pos_diff=1.0, max_cell_diff=1.0, **castep_args)¶

Bases: matscipy.socketcalc.QMClient

Subclass of Client for running CASTEP calculations.

Initial input files are written in .cell and .param formats, and subsequent communication is via sockets in REFTRAJ format.

extra_args(label=None)¶: Return list of additional command line arguments to be passed to client

postprocess(at, label)¶

Post-process results of calculation.

May be overrriden in subclasses to e.g. reverse sort order applied in preprocess().

preprocess(at, label, force_restart=False)¶

Prepare client for a calculation.

Starts client if this is the first task for it, or schedules a restart if new configuration is not compatible with the last one submitted to the queue (see is_compatible() method).

Many be extended in subclasses to e.g. sort the atoms by atomic number. If Atoms object needs to be changed, a copy should be returned rather than updating it inplace.

Returns (at, first_time).

write_input_files(at, label)¶

class matscipy.socketcalc.Client(client_id, exe, env=None, npj=1, ppn=1, block=None, corner=None, shape=None, jobname='socketcalc', rundir=None, fmt='REFTRAJ', parmode=None, mpirun='mpirun', mpirun_args=['-np'], logger=<matscipy.logger.Logger object>, max_pos_diff=1.0, max_cell_diff=0.001)¶

Bases: object

Represents a single Client job

Used by AtomsServer to start, restart and shutdown clients running on the Compute Nodes.

extra_args(label=None)¶: Return list of additional command line arguments to be passed to client

is_compatible(old_at, new_at, label)¶

Check if new_at and old_at are compatible.

Returns True if calculation can be continued, or False if client must be restarted before it can process new_at.

postprocess(at, label)¶

Post-process results of calculation.

May be overrriden in subclasses to e.g. reverse sort order applied in preprocess().

preprocess(at, label, force_restart=False)¶

Prepare client for a calculation.

Starts client if this is the first task for it, or schedules a restart if new configuration is not compatible with the last one submitted to the queue (see is_compatible() method).

Many be extended in subclasses to e.g. sort the atoms by atomic number. If Atoms object needs to be changed, a copy should be returned rather than updating it inplace.

Returns (at, first_time).

shutdown(block=True)¶

Request a client to shutdown.

If block=True, does not return until shutdown is complete. If block=False, waits for the client to shutdown in a new thread. Check self.waits_thread.isAlive() to see when shutdown has finished. (This function also returns a handle to the wait thread when block=False).

start(label=None)¶

Start an individual client.

Raises RuntimeError if this client is already running.

start_or_restart(at, label, restart=False)¶

Start or restart a client

If restart=True, wait for previous client to shutdown first. Calls write_input_files() followed by start().

wait_for_shutdown()¶

Block until a client has shutdown.

Typically called automatically by shutdown() or start_or_restart().

Shutdown should previously have been initiated by queuing a ‘shutdown’ or ‘restart’ request. Waits CLIENT_TIMEOUT for graceful shutdown. If client is still alive, a SIGTERM signal is sent. If this has had no effect after a further CLIENT_TIMEOUT, then a SIGKILL is sent. Does not return until the SIGKILL has taken effect.

This function also marks shutdown task as complete in servers’s input_q for this client.

write_input_files(at, label)¶

class matscipy.socketcalc.QMClient(client_id, exe, env=None, npj=1, ppn=1, block=None, corner=None, shape=None, jobname='socketcalc', rundir=None, fmt='REFTRAJ', parmode=None, mpirun='mpirun', mpirun_args=['-np'], logger=<matscipy.logger.Logger object>, max_pos_diff=1.0, max_cell_diff=0.001)¶

Bases: matscipy.socketcalc.Client

Abstract subclass of Client for QM calculations

is_compatible(old_at, new_at, label)¶

Check if new_at and old_at are compatible.

Returns True if calculation can be continued, or False if client must be restarted before it can process new_at.

class matscipy.socketcalc.QUIPClient(client_id, exe, env=None, npj=1, ppn=1, block=None, corner=None, shape=None, jobname='socketcalc', rundir=None, fmt='REFTRAJ', parmode=None, mpirun='mpirun', mpirun_args=['-np'], logger=<matscipy.logger.Logger object>, max_pos_diff=1.0, max_cell_diff=0.001, param_files=None)¶

Bases: matscipy.socketcalc.Client

Subclass of Client for running QUIP calculations.

Initial input files are written in extended XYZ format, and subsequent communication is via sockets, in either REFTRAJ or XYZ format.

write_input_files(at, label)¶

class matscipy.socketcalc.SocketCalculator(client, ip=None, atoms=None, port=0, logger=<matscipy.logger.Logger object>, bgq=False)¶

Bases: ase.calculators.calculator.Calculator

ASE-compatible calculator which communicates with remote force engines via sockets using a (synchronous) AtomsServer.

calculate(atoms, properties, system_changes)¶

Do the calculation.

properties: list of str: List of what needs to be calculated. Can be any combination of ‘energy’, ‘forces’, ‘stress’, ‘dipole’, ‘charges’, ‘magmom’ and ‘magmoms’.
system_changes: list of str: List of what has changed since last calculation. Can be any combination of these six: ‘positions’, ‘numbers’, ‘cell’, ‘pbc’, ‘initial_charges’ and ‘initial_magmoms’.

Subclasses need to implement this, but can ignore properties and system_changes if they want. Calculated properties should be inserted into results dictionary like shown in this dummy example:

self.results = {'energy': 0.0,
                'forces': np.zeros((len(atoms), 3)),
                'stress': np.zeros(6),
                'dipole': np.zeros(3),
                'charges': np.zeros(len(atoms)),
                'magmom': 0.0,
                'magmoms': np.zeros(len(atoms))}

The subclass implementation should first call this implementation to set the atoms attribute.

default_parameters = {}¶

implemented_properties = ['energy', 'forces', 'stress']¶

name = 'SocketCalculator'¶

shutdown()¶

class matscipy.socketcalc.VaspClient(client_id, exe, env=None, npj=1, ppn=1, block=None, corner=None, shape=None, jobname='socketcalc', rundir=None, fmt='REFTRAJ', parmode=None, mpirun='mpirun', mpirun_args=['-np'], logger=<matscipy.logger.Logger object>, max_pos_diff=1.0, max_cell_diff=0.001, **vasp_args)¶

Bases: matscipy.socketcalc.QMClient

Subclass of Client for running VASP calculations.

Initial input files are written in POSCAR, INCAR, POTCAR and KPOINTS formats, and subsequent communicatin is via sockets in REFTRAJ format.

postprocess(at, label)¶

Post-process results of calculation.

May be overrriden in subclasses to e.g. reverse sort order applied in preprocess().

preprocess(at, label, force_restart=False)¶

Prepare client for a calculation.

Starts client if this is the first task for it, or schedules a restart if new configuration is not compatible with the last one submitted to the queue (see is_compatible() method).

Many be extended in subclasses to e.g. sort the atoms by atomic number. If Atoms object needs to be changed, a copy should be returned rather than updating it inplace.

Returns (at, first_time).

write_input_files(at, label)¶

matscipy.socketcalc.pack_atoms_to_reftraj_str(at, label)¶

matscipy.socketcalc.pack_atoms_to_xyz_str(at, label)¶

matscipy.socketcalc.pack_results_to_reftraj_output_str(at)¶

matscipy.socketcalc.unpack_reftraj_output_str_to_results(data)¶

matscipy.socketcalc.unpack_reftraj_str_to_atoms(data)¶

matscipy.socketcalc.unpack_xyz_str_to_results(data)¶

matscipy.structure_identification module¶

matscipy.surface module¶

matscipy.surface.MillerDirection(v)¶: Special case of MillerIndex with type="direction" (the default)

class matscipy.surface.MillerIndex¶

Bases: numpy.ndarray

Representation of a three of four index Miller direction or plane

A MillerIndex can be constructed from vector or parsed from a string:

x = MillerIndex('-211')
y = MillerIndex('111', type='plane')
z = x.cross(y)
print x # prints "[-211]"
print y # prints "(111)", note round brackets denoting a plane
print z.latex()
assert(angle_between(x,y) == pi/2.)
assert(angle_between(y,z) == pi/2.)
assert(angle_between(x,z) == pi/2.)

all_brackets = ['[', ']', '<', '>', '(', ')', '{', '}']¶

angle(other)¶

as3()¶

as4()¶

brackets = {'direction': '[]', 'direction_family': '<>', 'plane': '()', 'plane_family': '{}'}¶

cosine(other)¶

cross(other)¶

hat()¶

latex()¶: Format this MillerIndex as a LaTeX string

norm()¶

normalised()¶

classmethod parse(s)¶

Parse a Miller index string

Negative indices can be denoted by:

leading minus sign, e.g. [11-2]
trailing b (for ‘bar’), e.g. 112b
LaTeX \bar{}, e.g. [11\bar{2}] (which renders as \([11\bar{2}]\) in LaTeX)

Leading or trailing brackets of various kinds are ignored. i.e. [001], {001}, (001), [001], <001>, 001 are all equivalent.

Returns an array of components (i,j,k) or (h,k,i,l)

plane_spacing(a)¶

simplified()¶

simplify()¶: Simplify by dividing through by greatest common denominator

matscipy.surface.MillerPlane(v)¶: Special case of MillerIndex with type="plane"

matscipy.surface.angle_between(a, b)¶: Angle between crystallographic directions between a=[ijk] and b=[lmn], in radians.

matscipy.surface.gcd(a, b)¶: Calculate the greatest common divisor of a and b

matscipy.surface.make_unit_slab(unit_cell, axes)¶

General purpose unit slab creation routine

Only tested with cubic unit cells.

Code translated from quippy.structures.unit_slab(): https://github.com/libAtoms/QUIP/blob/public/src/libAtoms/Structures.f95

Parameters

unit_cell (Atoms) – Atoms object containing primitive unit cell
axes (3x3 array) – Miller indices of desired slab, as columns

Returns

slab – Output slab, with axes aligned with x, y, z.

Return type

Atoms

Module contents¶

matscipy.has_parameter(name)¶

Test if a parameter has been provided in params.py.

Parameters: name (str) – Name of the parameter.
Returns: value – Returns True if parameter exists.
Return type: bool

matscipy.parameter(name, default=None, logger=<matscipy.logger.Logger object>)¶

Read parameter from params.py control file.

Parameters

name (str) – Name of the parameter.
default (optional) – Default value. Will be returned if parameter is not present.

Returns

Value of the parameter.

Return type

value

matscipy package¶

Subpackages¶

Submodules¶

matscipy.angle_distribution module¶

matscipy.atomic_strain module¶

matscipy.elasticity module¶

matscipy.hydrogenate module¶

matscipy.io module¶

matscipy.logger module¶

matscipy.neighbours module¶

matscipy.rings module¶

matscipy.socketcalc module¶

matscipy.structure_identification module¶

matscipy.surface module¶

Module contents¶