API Reference

System

Classes which specify the physical system.

Notes

This file is part of OpenFerro.

class openferro.system.System(lattice)

Bases: object

A class to define a physical system.

A system is a lattice with fields and a Hamiltonian. Fields are added to the system by the user. Interactions are added to the system by the user. Pointers to fields and interactions are stored in dictionaries.

lattice

The lattice of the system

Type:

BravaisLattice3D

__init__(lattice)

get_field_by_ID(ID)

Get a field by ID.

Parameters:

ID (str) – ID of the field

Returns:

The field with the given ID

Return type:

Field

Raises:

ValueError – If field with given ID does not exist

get_all_fields()

Get all fields in the system.

Returns:

All fields in the system

Return type:

list

get_all_SO3_fields()

get_all_non_SO3_fields()

move_fields_to_multi_devs(mesh: DeviceMesh)

Move all fields to given devices for parallelization.

Parameters:

mesh (DeviceMesh) – The device mesh to move the fields to

add_field(ID, ftype='scalar', dim=None, value=None, mass=1.0)

Add a predefined field to the system.

Parameters:

• ID (str) – ID of the field

• ftype (str, optional) – Type of the field. Can be ‘scalar’, ‘SO3’, ‘LocalStrain3D’, etc

• dim (int, optional) – Dimension of the field. Only used for Rn fields

• value (array-like, optional) – Initial value of the field. Will be broadcasted to the shape of the field

• mass (float or array-like, optional) – Mass of the field. When mass is a float, it will be broadcasted to the shape of the field

Returns:

The field with the given ID

Return type:

Field

Raises:

• ValueError – If field with this ID already exists or if ID is reserved

• ValueError – If field type is unknown

add_global_strain(value=None, mass=1)

Add a global strain to the system. Allow variable cell simulation.

Parameters:

• value (array-like, optional) – Initial value of the global strain

• mass (float, optional) – Effective mass of the global strain for the barostat

Returns:

The global strain field

Return type:

Field

Raises:

AssertionError – If value is provided but not a 6D vector

property interaction_dict

Get all interactions in the system.

Returns:

All interactions in the system

Return type:

dict

get_interaction_by_ID(interaction_ID)

Get an interaction by ID.

Parameters:

interaction_ID (str) – ID of the interaction

Returns:

The interaction with the given ID

Return type:

Interaction

Raises:

ValueError – If interaction with given ID does not exist

_add_interaction_sanity_check(ID)

Sanity check for adding an interaction.

Parameters:

ID (str) – ID of the interaction to check

Raises:

ValueError – If ID is reserved or already exists

add_dipole_dipole_interaction(ID, field_ID, prefactor=1.0, enable_jit=True)

Add the long-range dipole-dipole interaction term to the Hamiltonian.

Parameters:

• ID (str) – ID of the interaction

• field_ID (str) – ID of the field

• prefactor (float, optional) – Prefactor of the interaction

• enable_jit (bool, optional) – Whether to use JIT compilation

Returns:

The created interaction

Return type:

Interaction

add_dipole_onsite_interaction(ID, field_ID, K2, alpha, gamma, enable_jit=True)

add_dipole_interaction_1st_shell(ID, field_ID, j1, j2, enable_jit=True)

add_dipole_interaction_2nd_shell(ID, field_ID, j3, j4, j5, enable_jit=True)

add_dipole_interaction_3rd_shell(ID, field_ID, j6, j7, enable_jit=True)

add_dipole_efield_interaction(ID, field_ID, E, enable_jit=True)

add_homo_elastic_interaction(ID, field_ID, B11, B12, B44, enable_jit=True)

add_inhomo_elastic_interaction(ID, field_ID, B11, B12, B44, enable_jit=True)

add_pressure(pressure)

Add a pressure term (pV) to the Hamiltonian.

The ID of the interaction is reserved as ‘pV’. V is the volume of the system, which is calculated from the reference lattice vectors and the global strain.

Parameters:

pressure (float) – Pressure in bars

Returns:

The pV interaction

Return type:

Interaction

Raises:

ValueError – If pV term already exists or if gstrain field is invalid

add_homo_strain_dipole_interaction(ID, field_1_ID, field_2_ID, B1xx, B1yy, B4yz, enable_jit=True)

add_inhomo_strain_dipole_interaction(ID, field_1_ID, field_2_ID, B1xx, B1yy, B4yz, enable_jit=True)

add_cubic_anisotropy_interaction(ID, field_ID, K1, K2, enable_jit=True)

Add the cubic anisotropy interaction term.

Parameters:

• ID (str) – ID of the interaction

• field_ID (str) – ID of the field

• K1 (float) – First anisotropy constant

• K2 (float) – Second anisotropy constant

• enable_jit (bool, optional) – Whether to use JIT compilation

Returns:

The created interaction

Return type:

Interaction

_add_isotropic_exchange_interaction_by_rollers(ID, field_ID, coupling, rollers, enable_jit=True)

Add the isotropic exchange interaction term H=sum_{i~j} Jij*Si*Sj to the Hamiltonian.

Parameters:

• ID (str) – ID of the interaction

• field_ID (str) – ID of the field

• coupling (float) – Coupling constant

• rollers (list) – List of rolling functions for specifying the neighbouring relationship

• enable_jit (bool, optional) – Whether to use JIT compilation

Returns:

The created interaction

Return type:

Interaction

add_isotropic_exchange_interaction_1st_shell(ID, field_ID, coupling, enable_jit=True)

Add the first shell isotropic exchange interaction term.

The first shell is defined in lattice class.

Parameters:

• ID (str) – ID of the interaction

• field_ID (str) – ID of the field

• coupling (float) – Coupling constant

• enable_jit (bool, optional) – Whether to use JIT compilation

Returns:

The created interaction

Return type:

Interaction

add_isotropic_exchange_interaction_2nd_shell(ID, field_ID, coupling, enable_jit=True)

Add the second shell isotropic exchange interaction term.

The second shell is defined in lattice class.

Parameters:

• ID (str) – ID of the interaction

• field_ID (str) – ID of the field

• coupling (float) – Coupling constant

• enable_jit (bool, optional) – Whether to use JIT compilation

Returns:

The created interaction

Return type:

Interaction

add_isotropic_exchange_interaction_3rd_shell(ID, field_ID, coupling, enable_jit=True)

Add the third shell isotropic exchange interaction term.

The third shell is defined in lattice class.

Parameters:

• ID (str) – ID of the interaction

• field_ID (str) – ID of the field

• coupling (float) – Coupling constant

• enable_jit (bool, optional) – Whether to use JIT compilation

Returns:

The created interaction

Return type:

Interaction

add_isotropic_exchange_interaction_4th_shell(ID, field_ID, coupling, enable_jit=True)

Add the fourth shell isotropic exchange interaction term.

The fourth shell is defined in lattice class.

Parameters:

• ID (str) – ID of the interaction

• field_ID (str) – ID of the field

• coupling (float) – Coupling constant

• enable_jit (bool, optional) – Whether to use JIT compilation

Returns:

The created interaction

Return type:

Interaction

add_self_interaction(ID, field_ID, energy_engine, parameters=None, enable_jit=True)

Add a custom self-interaction term to the Hamiltonian.

Parameters:

• ID (str) – ID of the interaction

• field_ID (str) – ID of the field

• energy_engine (callable) – A function that takes the field as input and returns the interaction energy

• parameters (array-like, optional) – Parameters for the interaction

• enable_jit (bool, optional) – Whether to use JIT compilation

Returns:

The created interaction

Return type:

Interaction

add_mutual_interaction(ID, field_1_ID, field_2_ID, energy_engine, parameters=None, enable_jit=True)

Add a custom mutual interaction term to the Hamiltonian.

Parameters:

• ID (str) – ID of the interaction

• field_1_ID (str) – ID of the first field

• field_2_ID (str) – ID of the second field

• energy_engine (callable) – A function that takes the fields as input and returns the interaction energy

• parameters (array-like, optional) – Parameters for the interaction

• enable_jit (bool, optional) – Whether to use JIT compilation

Returns:

The created interaction

Return type:

Interaction

add_triple_interaction(ID, field_1_ID, field_2_ID, field_3_ID, energy_engine, parameters=None, enable_jit=True)

Add a custom triple interaction term to the Hamiltonian.

Parameters:

• ID (str) – ID of the interaction

• field_1_ID (str) – ID of the first field

• field_2_ID (str) – ID of the second field

• field_3_ID (str) – ID of the third field

• energy_engine (callable) – A function that takes the fields as input and returns the interaction energy

• parameters (array-like, optional) – Parameters for the interaction

• enable_jit (bool, optional) – Whether to use JIT compilation

Returns:

The created interaction

Return type:

Interaction

calc_energy_by_ID(interaction_ID)

Calculate the energy of an interaction by ID.

Parameters:

interaction_ID (str) – ID of the interaction

Returns:

The energy of the interaction

Return type:

float

Raises:

ValueError – If interaction with given ID does not exist

calc_force_by_ID(interaction_ID)

Calculate the gradient force from an interaction by ID.

Parameters:

interaction_ID (str) – ID of the interaction

Returns:

The gradient force of the interaction

Return type:

array

Raises:

ValueError – If interaction with given ID does not exist

calc_total_self_energy()

Calculate the total self-interaction energy.

Returns:

Total self-interaction energy

Return type:

float

calc_total_mutual_interaction()

Calculate the total mutual interaction energy.

Returns:

Total mutual interaction energy

Return type:

float

calc_total_triple_interaction()

Calculate the total triple interaction energy.

Returns:

Total triple interaction energy

Return type:

float

calc_total_potential_energy()

Calculate the total potential energy of the system.

Calculate the total potential energy of the system. Total potential energy is the sum of self-interaction energy, mutual interaction energy, and triple interaction energy.

calc_total_kinetic_energy()

Calculate the total kinetic energy of the system.

calc_temp_by_ID(ID)

Calculate the temperature of a field.

calc_excess_stress()

Get instantaneous stress - applied stress (e.g. from hydrostatic pressure)

update_force_from_self_interaction(profile=False)

update the gradient force felt by each field from self interactions.

update_force_from_mutual_interaction(profile=False)

update the gradient force felt by each field from mutual interactions.

update_force_from_triple_interaction(profile=False)

update the gradient force felt by each field from triple interactions.

update_force(profile=False)

update the gradient force felt by each field from all interactions.

class openferro.system.RingPolymerSystem(lattice, nbeads=1)

Bases: System

A class to define a ring polymer system for path-integral molecular dynamics simulations.

__init__(lattice, nbeads=1)

Fields

Classes which define the fields on the lattice.

class openferro.field.Field(lattice, ID: str)

Bases: object

Template class to define a field on a lattice.

__init__(lattice, ID: str)

Initialize a field.

Parameters:

• lattice (BravaisLattice3D) – Lattice object

• ID (str) – ID of the field

set_values(values)

get_values()

Get the values of the field.

Returns:

Values of the field

Return type:

array_like

Raises:

ValueError – If field has no values set

property size

set_mass(mass)

get_mass()

set_velocity(velocity)

get_velocity()

init_velocity(mode='zero', temperature=None)

set_force(force)

get_force()

zero_force()

accumulate_force(force)

get_kinetic_energy()

get_temperature()

compare_shape(x, y)

compare_sharding(x, y)

to_multi_devs(mesh: DeviceMesh)

set_integrator(integrator_class, dt, **kwargs)

Set the integrator according to the given integrator class. Set the time step. To be implemented by the subclasses.

Parameters:

• integrator_class (str) – Class of integrator to use

• dt (float) – Time step

• **kwargs – Additional arguments passed to integrator

set_custom_integrator(integrator)

class openferro.field.FieldRn(lattice, ID, dim, unit=None)

Bases: Field

R^n field on a lattice. Values are stored as n-dimensional vectors.

__init__(lattice, ID, dim, unit=None)

Initialize a field.

Parameters:

• lattice (BravaisLattice3D) – Lattice object

• ID (str) – ID of the field

property mean

Calculate the average of the field over the lattice.

Returns:

Mean value of field

Return type:

array_like

property var

Calculate the variance of the field over the lattice.

Returns:

Variance of field

Return type:

array_like

set_local_value(loc, value)

Set the value of the field at a given location.

Parameters:

• loc (tuple) – Location tuple with length equal to lattice dimension

• value (array_like) – Value to set at location

set_integrator(integrator_class, dt, temp=None, tau=None)

Set the integrator according to the given integrator class. Set the time step.

Parameters:

• integrator_class (str) – Integrator class

• dt (float) – Time step

• temp (float, optional) – Temperature for isothermal integrator

• tau (float, optional) – Relaxation time for Langevin integrator

Raises:

ValueError – If integrator class not supported or required parameters missing

class openferro.field.FieldScalar(lattice, ID, unit=None)

Bases: FieldRn

Scalar field. Values are stored as scalars. Example: mass, density, any onsite constant, etc.

__init__(lattice, ID, unit=None)

Initialize a field.

Parameters:

• lattice (BravaisLattice3D) – Lattice object

• ID (str) – ID of the field

class openferro.field.FieldR3(lattice, ID, unit=None)

Bases: FieldRn

3D vector field. Values are stored as 3-dimensional vectors. Example: flexible dipole moment.

__init__(lattice, ID, unit=None)

Initialize a field.

Parameters:

• lattice (BravaisLattice3D) – Lattice object

• ID (str) – ID of the field

class openferro.field.FieldSO3(lattice, ID, unit=None)

Bases: FieldRn

3D vector field with fixed magnitude and flexible orientation. Values are stored as 3-dimensional vectors. Example: rigid atomistic spin, molecular orientation, etc.

__init__(lattice, ID, unit=None)

Initialize a field.

Parameters:

• lattice (BravaisLattice3D) – Lattice object

• ID (str) – ID of the field

set_magnitude(magnitude)

get_magnitude()

perturb(sigma)

normalize()

init_velocity(mode='zero', temperature=None)

to_multi_devs(mesh: DeviceMesh)

set_integrator(integrator_class, dt, temp=None, alpha=None)

Set the integrator according to the given integrator class. Set the time step.

Parameters:

• integrator_class (str) – Integrator class

• dt (float) – Time step

• temp (float, optional) – Temperature for isothermal integrator

• alpha (float, optional) – Gilbert damping constant for Landau-Lifshitz equation of motion

Raises:

ValueError – If integrator class not supported or required parameters missing

class openferro.field.LocalStrain3D(lattice, ID)

Bases: FieldRn

Strain field on 3D lattice are separated into local contribution (local strain field) and global contribution (homogeneous strain associated to the supercell). The local strain field is encoded by the local displacement vector v_i(R)/a_i (a_i: the lattice vector) associated with each lattice site at R.

__init__(lattice, ID)

Initialize a field.

Parameters:

• lattice (BravaisLattice3D) – Lattice object

• ID (str) – ID of the field

static get_local_strain_symmetric(values)

Calculate the local strain field from the local displacement field.

Parameters:

values (array_like) – Local displacement field values

Returns:

Local strain field with shape (l1, l2, l3, 6)

Return type:

array_like

static get_local_strain(values)

Calculate the local strain field from the local displacement field. Implemented according to Physical Review B 52.9 (1995): 6301.

Parameters:

values (array_like) – Local displacement field values

Returns:

Local strain field with shape (l1, l2, l3, 6)

Return type:

array_like

class openferro.field.GlobalStrain(lattice, ID)

Bases: Field

The homogeneous strain is represented by the strain tensor with Voigt convention, which is a 6-dimensional vector.

__init__(lattice, ID)

Initialize a field.

Parameters:

• lattice (BravaisLattice3D) – Lattice object

• ID (str) – ID of the field

to_multi_devs(mesh: DeviceMesh)

get_excess_stress()

set_integrator(integrator_class, dt, temp=None, tau=None, freeze_x=False, freeze_y=False, freeze_z=False)

Set the integrator according to the given integrator class. Set the time step.

Parameters:

• integrator_class (str) – Integrator class

• dt (float) – Time step

• temp (float, optional) – Temperature for isothermal integrator

• tau (float, optional) – Relaxation time for Langevin integrator

• freeze_x (bool, optional) – Whether to freeze x-component of strain, by default False

• freeze_y (bool, optional) – Whether to freeze y-component of strain, by default False

• freeze_z (bool, optional) – Whether to freeze z-component of strain, by default False

Raises:

ValueError – If integrator class not supported or required parameters missing

Lattice

Classes to define the Bravais lattices

Notes

This file is part of OpenFerro.

class openferro.lattice.BravaisLattice3D(l1, l2, l3, a1=None, a2=None, a3=None, pbc=True)

Bases: object

Base class to represent a 3D Bravais lattice

Parameters:

• l1 (int) – Size in first dimension

• l2 (int) – Size in second dimension

• l3 (int) – Size in third dimension

• a1 (array_like, optional) – First primitive vector, defaults to [1,0,0]

• a2 (array_like, optional) – Second primitive vector, defaults to [0,1,0]

• a3 (array_like, optional) – Third primitive vector, defaults to [0,0,1]

• pbc (bool, optional) – Whether to use periodic boundary conditions, defaults to True

__init__(l1, l2, l3, a1=None, a2=None, a3=None, pbc=True)

property nsites

property unit_volume

Returns the volume of the unit cell

Returns:

Volume of unit cell

Return type:

float

property ref_volume

property latt_vec

Returns lattice vectors.

Returns:

3x3 array of lattice vectors

Return type:

ndarray

property reciprocal_latt_vec

Calculates reciprocal lattice vectors.

Returns:

3x3 array of reciprocal lattice vectors

Return type:

ndarray

class openferro.lattice.SimpleCubic3D(l1, l2, l3, a1=None, a2=None, a3=None, pbc=True)

Bases: BravaisLattice3D

A class to represent a simple cubic lattice

Parameters:

• l1 (int) – Size in first dimension

• l2 (int) – Size in second dimension

• l3 (int) – Size in third dimension

• a1 (array_like, optional) – First primitive vector

• a2 (array_like, optional) – Second primitive vector

• a3 (array_like, optional) – Third primitive vector

• pbc (bool, optional) – Whether to use periodic boundary conditions, defaults to True

Notes

Coordination number is 6. First shell has 6 neighbours. Second shell has 12 neighbours. Third shell has 8 neighbours.

__init__(l1, l2, l3, a1=None, a2=None, a3=None, pbc=True)

_1st_shell_roller()

Return a list of rolling functions for moving a site to all sites in the first shell of the lattice.

Returns:

List of rolling functions. The number of rolling functions is half the shell coordination number. Used to calculate interactions between a site and its neighbours.

Return type:

list

_2nd_shell_roller()

Return a list of rolling functions for moving a site to all sites in the second shell of the lattice.

Returns:

List of rolling functions. The number of rolling functions is half the shell coordination number. Used to calculate interactions between a site and its neighbours.

Return type:

list

_3rd_shell_roller()

Return a list of rolling functions for moving a site to all sites in the third shell of the lattice.

Returns:

List of rolling functions. The number of rolling functions is half the shell coordination number. Used to calculate interactions between a site and its neighbours.

Return type:

list

class openferro.lattice.BodyCenteredCubic3D(l1, l2, l3, a1=None, a2=None, a3=None, pbc=True)

Bases: SimpleCubic3D

A class to represent a body-centered cubic lattice

Parameters:

• l1 (int) – Size in first dimension

• l2 (int) – Size in second dimension

• l3 (int) – Size in third dimension

• a1 (array_like, optional) – First primitive vector

• a2 (array_like, optional) – Second primitive vector

• a3 (array_like, optional) – Third primitive vector

• pbc (bool, optional) – Whether to use periodic boundary conditions, defaults to True

Notes

Coordination number is 8. First shell has 8 neighbours. Second shell has 6 neighbours. Third shell has 12 neighbours. Fourth shell has 24 neighbours.

__init__(l1, l2, l3, a1=None, a2=None, a3=None, pbc=True)

_1st_shell_roller()

Return a list of rolling functions for moving a site to all sites in the first shell of the lattice.

Returns:

List of rolling functions. The number of rolling functions is half the shell coordination number. Used to calculate interactions between a site and its neighbours.

Return type:

list

_2nd_shell_roller()

Return a list of rolling functions for moving a site to all sites in the second shell of the lattice.

Returns:

List of rolling functions. The number of rolling functions is half the shell coordination number. Used to calculate interactions between a site and its neighbours.

Return type:

list

_3rd_shell_roller()

Return a list of rolling functions for moving a site to all sites in the third shell of the lattice.

Returns:

List of rolling functions. The number of rolling functions is half the shell coordination number. Used to calculate interactions between a site and its neighbours.

Return type:

list

_4th_shell_roller()

Return a list of rolling functions for moving a site to all sites in the fourth shell of the lattice.

Returns:

List of rolling functions. The number of rolling functions is half the shell coordination number. Used to calculate interactions between a site and its neighbours.

Return type:

list

class openferro.lattice.FaceCenteredCubic3D(l1, l2, l3, a1=None, a2=None, a3=None, pbc=True)

Bases: BodyCenteredCubic3D

A class to represent a face-centered cubic lattice

Parameters:

• l1 (int) – Size in first dimension

• l2 (int) – Size in second dimension

• l3 (int) – Size in third dimension

• a1 (array_like, optional) – First primitive vector

• a2 (array_like, optional) – Second primitive vector

• a3 (array_like, optional) – Third primitive vector

• pbc (bool, optional) – Whether to use periodic boundary conditions, defaults to True

Notes

Coordination number is 12. First shell has 12 neighbours. Second shell has 6 neighbours. Third shell has 24 neighbours.

__init__(l1, l2, l3, a1=None, a2=None, a3=None, pbc=True)

class openferro.lattice.Hexagonal3D(l1, l2, l3, a1=None, a2=None, a3=None, pbc=True)

Bases: BravaisLattice3D

A class to represent a hexagonal lattice

Parameters:

• l1 (int) – Size in first dimension

• l2 (int) – Size in second dimension

• l3 (int) – Size in third dimension

• a1 (array_like, optional) – First primitive vector

• a2 (array_like, optional) – Second primitive vector

• a3 (array_like, optional) – Third primitive vector

• pbc (bool, optional) – Whether to use periodic boundary conditions, defaults to True

__init__(l1, l2, l3, a1=None, a2=None, a3=None, pbc=True)

class openferro.lattice.BravaisLattice2D(l1, l2, a1=None, a2=None, pbc=True)

Bases: object

Base class to represent a 2D Bravais lattice

Parameters:

• l1 (int) – Size in first dimension

• l2 (int) – Size in second dimension

• a1 (array_like, optional) – First primitive vector

• a2 (array_like, optional) – Second primitive vector

• pbc (bool, optional) – Whether to use periodic boundary conditions, defaults to True

__init__(l1, l2, a1=None, a2=None, pbc=True)

property unit_area

Returns the area of the unit cell

Returns:

Area of unit cell

Return type:

float

property ref_area

property latt_vec

Returns lattice vectors.

Returns:

2x2 array of lattice vectors

Return type:

ndarray

property reciprocal_latt_vec

Calculates reciprocal lattice vectors.

Returns:

2x2 array of reciprocal lattice vectors

Return type:

ndarray

class openferro.lattice.SimpleSquare2D(l1, l2, a1=None, a2=None, pbc=True)

Bases: BravaisLattice2D

A class to represent a simple square lattice

Parameters:

• l1 (int) – Size in first dimension

• l2 (int) – Size in second dimension

• a1 (array_like, optional) – First primitive vector

• a2 (array_like, optional) – Second primitive vector

• pbc (bool, optional) – Whether to use periodic boundary conditions, defaults to True

Notes

Coordination number is 4. First shell has 4 neighbours. Second shell has 4 neighbours. Third shell has 4 neighbours.

__init__(l1, l2, a1=None, a2=None, pbc=True)

_1st_shell_roller()

Return a list of rolling functions for moving a site to all sites in the first shell of the lattice.

Returns:

List of rolling functions. The number of rolling functions is half the shell coordination number. Used to calculate interactions between a site and its neighbours.

Return type:

list

_2nd_shell_roller()

Return a list of rolling functions for moving a site to all sites in the second shell of the lattice.

Returns:

List of rolling functions. The number of rolling functions is half the shell coordination number. Used to calculate interactions between a site and its neighbours.

Return type:

list

_3rd_shell_roller()

Return a list of rolling functions for moving a site to all sites in the third shell of the lattice.

Returns:

List of rolling functions. The number of rolling functions is half the shell coordination number. Used to calculate interactions between a site and its neighbours.

Return type:

list

class openferro.lattice.Hexagonal2D(l1, l2, a1=None, a2=None)

Bases: SimpleSquare2D

A class to represent a hexagonal lattice

Parameters:

• l1 (int) – Size in first dimension

• l2 (int) – Size in second dimension

• a1 (array_like, optional) – First primitive vector

• a2 (array_like, optional) – Second primitive vector

Notes

Coordination number is 6.

__init__(l1, l2, a1=None, a2=None)

Hamiltonian

Elastic

Functions that define a term in the elastic energy.

These functions will be added into <class interaction> for automatic differentiation.

Notes

This file is part of OpenFerro.

openferro.engine.elastic.homo_elastic_energy(global_strain, parameters)

Returns the homogeneous elastic energy of a strain field.

Parameters:

• global_strain (jnp.ndarray) – The global strain of a supercell, shape=(6)

• parameters (jnp.ndarray) – The parameters of the energy function

Returns:

The homogeneous elastic energy

Return type:

jnp.ndarray

openferro.engine.elastic.pV_energy(global_strain, parameters)

Returns pressure * (volume - reference volume)

Parameters:

• global_strain (jnp.ndarray) – The global strain of a supercell, shape=(6)

• parameters (jnp.ndarray) – The parameters of the energy function

Returns:

The pV energy

Return type:

jnp.ndarray

openferro.engine.elastic.inhomo_elastic_energy(local_displacement, parameters)

Returns the inhomogeneous elastic energy of a strain field.

Parameters:

• local_displacement (jnp.ndarray) – The local displacement field, shape=(nx, ny, nz, 3)

• parameters (tuple) – The parameters of the energy function containing: - B11 (float): elastic constant B11 - B12 (float): elastic constant B12 - B44 (float): elastic constant B44

Returns:

The total elastic energy (homogeneous + inhomogeneous)

Return type:

jnp.ndarray

Ferroelectric

Functions that define a ferroelectric term in the Hamiltonian. They will be added into <class interaction> for automatic differentiation.

openferro.engine.ferroelectric.self_energy_onsite_isotropic(field, parameters)

Returns the isotropic self-energy of a 3D field. See Eq.(2-3) in [Zhong, W., David Vanderbilt, and K. M. Rabe. Physical Review B 52.9 (1995): 6301.] for meaning of the parameters.

Parameters:

• field (jnp.array) – The field to calculate the energy

• parameters (jax.numpy array) – The parameters of the energy function

Returns:

The energy of the field

Return type:

jnp.array

openferro.engine.ferroelectric.self_energy_onsite_scalar(field, parameters)

Returns the self-energy of a scalar field. E= sum_i E_i. (sum over the lattice sites i) E_i = k_2 * u_i^2 + alpha * u_i^4

Parameters:

• field (jnp.array) – The field to calculate the energy

• parameters (jax.numpy array) – The parameters of the energy function

Returns:

The energy of the field

Return type:

jnp.array

openferro.engine.ferroelectric.short_range_1stnn_isotropic_scalar(field, parameters)

Returns the short-range interaction of nearest neighbors for a R^3 field defined on a lattice with periodic boundary conditions.

Parameters:

• field (jnp.array) – The field to calculate the energy

• parameters (jax.numpy array) – The parameters of the energy function

Returns:

The energy of the field

Return type:

jnp.array

openferro.engine.ferroelectric.short_range_1stnn_isotropic(field, parameters)

Returns the short-range interaction of nearest neighbors for a R^3 field defined on a isotropic lattice with periodic boundary conditions.

Parameters:

• field (jnp.array) – The field to calculate the energy

• parameters (jax.numpy array) – The parameters of the energy function

Returns:

The energy of the field

Return type:

jnp.array

openferro.engine.ferroelectric.short_range_1stnn_anisotropic(field, parameters)

Returns the short-range interaction of nearest neighbors for a R^3 field defined on a anisotropic lattice with periodic boundary conditions.

Parameters:

• field (jnp.array) – The field to calculate the energy

• parameters (jax.numpy array) – The parameters of the energy function

Returns:

The energy of the field

Return type:

jnp.array

openferro.engine.ferroelectric.short_range_2ednn_isotropic(field, parameters)

Returns the short-range interaction of nearest neighbors for a R^3 field defined on a lattice with periodic boundary conditions.

Parameters:

• field (jnp.array) – The field to calculate the energy

• parameters (jax.numpy array) – The parameters of the energy function

Returns:

The energy of the field

Return type:

jnp.array

openferro.engine.ferroelectric.get_short_range_3rdnn_isotropic()

Returns the engine of short-range interaction of third nearest neighbors for a R^3 field defined on a lattice with periodic boundary conditions.

Returns:

The interaction function

Return type:

function

openferro.engine.ferroelectric.homo_strain_dipole_interaction(global_strain, dipole_field, parameters)

Returns the homogeneous strain dipole interaction energy.

Parameters:

• global_strain (jnp.array) – Shape=(6), the global strain of a supercell

• dipole_field (jnp.array) – Shape=(nx, ny, nz, 3), the dipole field

• parameters (jax.numpy array) – The parameters of the energy function containing B1xx : Elastic constant B1xx B1yy : Elastic constant B1yy B4yz : Elastic constant B4yz

Returns:

The homogeneous strain dipole interaction energy

Return type:

jnp.array

openferro.engine.ferroelectric.get_inhomo_strain_dipole_interaction(enable_jit=True)

Returns the inhomogeneous strain dipole interaction function.

Parameters:

enable_jit (bool, optional) – Whether to enable JIT compilation, by default True

Returns:

The interaction function

Return type:

function

openferro.engine.ferroelectric.dipole_efield_interaction(field, parameters)

Returns the dipole-electric field interaction energy.

Parameters:

• field (jnp.array) – The field to calculate the energy

• parameters (jax.numpy array) – The parameters of the energy function

Returns:

The interaction energy

Return type:

jnp.array

Magnetic

Functions that define a term in the magnetic Hamiltonian. They will be added into <class interaction> for automatic differentiation.

openferro.engine.magnetic.get_isotropic_exchange_energy_engine(rollers)

Returns the exchange energy engine for a R^3 field defined on a lattice with periodic boundary conditions.

Parameters:

rollers (list) – List of jnp.roll functions specifying the neighbors

Returns:

Energy engine function

Return type:

callable

openferro.engine.magnetic.cubic_anisotropy_energy(field, parameters)

Returns the anisotropy energy of the field.

Parameters:

• field (ndarray) – The magnetic field

• parameters (ndarray) – Array containing K1 and K2 anisotropy constants

Returns:

The anisotropy energy: -K1*(mx^2*my^2 + my^2*mz^2 + mx^2*mz^2) - K2*mx^2*my^2*mz^2

Return type:

float

openferro.engine.magnetic.Dzyaloshinskii_Moriya_energy(field, parameters)

Returns the Dzyaloshinskii-Moriya energy of the field.

Parameters:

• field (ndarray) – The magnetic field

• parameters (ndarray) – Array of parameters

Returns:

The Dzyaloshinskii-Moriya energy

Return type:

float

openferro.engine.magnetic.external_field_energy(field, parameters)

Returns the external field energy of the field.

Parameters:

• field (ndarray) – The magnetic field

• parameters (ndarray) – Array containing the external field B_ext

Returns:

The external field energy: -field·B_ext

Return type:

float

Ewald

Functions for Ewald summation

openferro.engine.ewald.get_dipole_dipole_ewald(latt, sharding=None)

Returns the function to calculate the energy of dipole-dipole interaction.

Implemented according to Sec.5.3 of “Wang, D., et al. ‘Ewald summation for ferroelectric perovksites with charges and dipoles.’ Computational Materials Science 162 (2019): 314-321.”

Parameters:

• latt (Lattice) – The lattice object containing size and lattice vectors

• sharding (jax.sharding.Sharding, optional) – Sharding specification for distributed arrays

Returns:

Function that calculates dipole-dipole interaction energy

Return type:

callable

openferro.engine.ewald.dipole_dipole_ewald_plain(field, parameters)

Brute-force Ewald summation for dipole-dipole interaction.

For benchmarking purpose only.

Parameters:

• field (ndarray) – The field values, shape=(l1, l2, l3, 3)

• parameters (dict) –

  Dictionary containing:

  a1float

  First lattice vector

  a2float

  Second lattice vector

  a3float

  Third lattice vector

  Z_starfloat

  Born effective charge

  epsilon_inffloat

  High-frequency dielectric constant

Returns:

The dipole-dipole interaction energy

Return type:

float

Integrators

Base

Base class for integrators.

This file is part of OpenFerro.

class openferro.integrator.base.Integrator(dt)

Bases: object

Base class for integrators.

Parameters:

dt (float) – Time step size

__init__(dt)

step(field, force_updater=None)

Update the field by one time step.

In most cases, the force will be updated for all fields in one setting, before any integrator is called. So the force_updater is not necessary in most cases. However, for some implicit integrators, the force_updater is needed.

Parameters:

• field (Field) – The field to be updated

• force_updater (callable, optional) – A function that updates the force of fields

Returns:

The updated field

Return type:

Field

Molecular Dynamics

Integrators for unconstrained molecular dynamics.

This file is part of OpenFerro.

class openferro.integrator.md.GradientDescentIntegrator(dt)

Bases: Integrator

Gradient descent integrator.

Parameters:

dt (float) – Time step size

__init__(dt)

step(field)

Update the field by one time step.

Parameters:

field (Field) – The field to be updated

Returns:

The updated field

Return type:

Field

class openferro.integrator.md.GradientDescentIntegrator_Strain(dt, freeze_x=False, freeze_y=False, freeze_z=False)

Bases: GradientDescentIntegrator

Gradient descent integrator for global strain.

Parameters:

• dt (float) – Time step size

• freeze_x (bool, optional) – Whether to freeze motion in x direction

• freeze_y (bool, optional) – Whether to freeze motion in y direction

• freeze_z (bool, optional) – Whether to freeze motion in z direction

__init__(dt, freeze_x=False, freeze_y=False, freeze_z=False)

class openferro.integrator.md.LeapFrogIntegrator(dt)

Bases: Integrator

Leapfrog integrator.

Parameters:

dt (float) – Time step size

__init__(dt)

step(field)

Update the field by one time step.

Parameters:

field (Field) – The field to be updated

Returns:

The updated field

Return type:

Field

class openferro.integrator.md.LeapFrogIntegrator_Strain(dt, freeze_x=False, freeze_y=False, freeze_z=False)

Bases: LeapFrogIntegrator

Leapfrog integrator for global strain.

Parameters:

• dt (float) – Time step size

• freeze_x (bool, optional) – Whether to freeze motion in x direction

• freeze_y (bool, optional) – Whether to freeze motion in y direction

• freeze_z (bool, optional) – Whether to freeze motion in z direction

__init__(dt, freeze_x=False, freeze_y=False, freeze_z=False)

class openferro.integrator.md.LangevinIntegrator(dt, temp, tau)

Bases: Integrator

Langevin integrator as in J. Phys. Chem. A 2019, 123, 28, 6056-6079. ABOBA scheme: exp(i L dt) = exp(i Lx dt/2)exp(i Lt dt)exp(i Lx dt/2)exp(i Lp dt)

Lx/2: half-step position update (_step_x) Lt: velocity update from damping and noise (_step_t) Lp: full-step velocity update (_step_p)

Parameters:

• dt (float) – Time step size

• temp (float) – Temperature

• tau (float) – Relaxation time

__init__(dt, temp, tau)

get_noise(key, field)

Generate random noise for the Langevin dynamics.

Parameters:

• key (jax.random.PRNGKey) – Random number generator key

• field (Field) – The field to generate noise for

Returns:

Random noise array

Return type:

jax.Array

step(key, field)

Update the field by one time step.

Parameters:

• key (jax.random.PRNGKey) – Random number generator key

• field (Field) – The field to be updated

Returns:

The updated field

Return type:

Field

class openferro.integrator.md.LangevinIntegrator_Strain(dt, temp, tau, freeze_x=False, freeze_y=False, freeze_z=False)

Bases: LangevinIntegrator

Langevin integrator for global strain.

Parameters:

• dt (float) – Time step size

• temp (float) – Temperature

• tau (float) – Relaxation time

• freeze_x (bool, optional) – Whether to freeze motion in x direction

• freeze_y (bool, optional) – Whether to freeze motion in y direction

• freeze_z (bool, optional) – Whether to freeze motion in z direction

__init__(dt, temp, tau, freeze_x=False, freeze_y=False, freeze_z=False)

class openferro.integrator.md.OverdampedLangevinIntegrator(dt, temp, tau)

Bases: Integrator

Overdamped Langevin integrator.

Parameters:

• dt (float) – Time step size

• temp (float) – Temperature

• tau (float) – Relaxation time

__init__(dt, temp, tau)

Landau-Lifshitz-Gilbert

Integrators for Landau-Lifshitz equations of motion. Fields are constrained to have fixed magnitude.

This file is part of OpenFerro.

class openferro.integrator.llg.ConservativeLLIntegrator(dt, gamma=None)

Bases: Integrator

[For testing purposes: naive geometric integrator that may cause increasement in energy.] Adiabatic spin precession, i.e. Landau-Lifshitz equation of motion without dissipative damping term.

The equation of motion is: dM/dt = -gamma M x B

Parameters:

• dt (float) – Time step size

• gamma (float, optional) – Gyromagnetic ratio. If None, uses electron gyromagnetic ratio

_step_x(M, B, gamma, dt)

Update the orientation of the magnetization M by rotating it around the magnetic field B

Parameters:

• M (jax.Array) – The magnetization (shape=(*, 3))

• B (jax.Array) – The magnetic field (shape=(*, 3))

• gamma (float) – The gyromagnetic ratio

• dt (float) – The time step

Returns:

The updated magnetization (shape=(*, 3))

Return type:

jax.Array

__init__(dt, gamma=None)

step(field, force_updater=None)

Update the field by one time step.

Parameters:

• field (Field) – The field to be updated

• force_updater (callable, optional) – A function that updates the force of fields

Returns:

The updated field

Return type:

Field

class openferro.integrator.llg.LLIntegrator(dt, alpha, gamma=None)

Bases: ConservativeLLIntegrator

[For testing purposes: naive geometric integrator.] Landau-Lifshitz equation of motion.

See Eriksson, Olle, et al. Atomistic spin dynamics: foundations and applications. Oxford university press, 2017, Sec.7.4.5 for details.

The equation of motion is

dM/dt = -gamma M x B - (gamma * alpha / |M|) * M x (M x B)

Here gamma=(gyromagnetic ratio)/ (1+alpha^2) is the renormalized gyromagnetic ratio for simulating LLG equation in Landau-Lifshitz form.

Parameters:

• dt (float) – Time step size

• alpha (float) – Gilbert damping constant

• gamma (float, optional) – Gyromagnetic ratio. If None, uses electron gyromagnetic ratio

_step_b(M, B, alpha, Ms)

Update the magnetic field B by adding the term ( alpha / |M|) * (M x B)

Parameters:

• M (jax.Array) – The magnetization (shape=(*, 3))

• B (jax.Array) – The magnetic field (shape=(*, 3))

• alpha (float) – The Gilbert damping constant

• Ms (jax.Array) – The magnitude of the magnetization (shape=(*,))

Returns:

The updated magnetic field (shape=(*, 3))

Return type:

jax.Array

__init__(dt, alpha, gamma=None)

step(field, force_updater=None)

Update the field by one time step.

Parameters:

• field (Field) – The field to be updated

• force_updater (callable, optional) – A function that updates the force of fields

Returns:

The updated field

Return type:

Field

class openferro.integrator.llg.LLLangevinIntegrator(dt, temp, alpha, gamma=None)

Bases: LLIntegrator

[For testing purposes: naive geometric integrator.] Stochastic Landau-Lifshitz equation of motion.

See Eriksson, Olle, et al. Atomistic spin dynamics: foundations and applications. Oxford university press, 2017, Sec.7.4.5 for details.

The equation of motion is: dM/dt = -gamma M x (B + b) - (gamma * alpha / |M|) * M x (M x (B + b))

b is the stochastic force <b_i_alpha(t) b_j_beta(s)> = 2 * D * delta_ij * delta_alpha_beta * delta(t-s) A steady Boltzmann state requires D= alpha/(1+alpha^2) * kbT / gamma / |m|

Parameters:

• dt (float) – Time step size

• temp (float) – Temperature

• alpha (float) – Gilbert damping constant

• gamma (float, optional) – Gyromagnetic ratio. If None, uses electron gyromagnetic ratio

__init__(dt, temp, alpha, gamma=None)

get_noise(key, field)

Generate random noise for the Langevin dynamics.

Parameters:

• key (jax.random.PRNGKey) – Random number generator key

• field (Field) – The field to generate noise for

Returns:

Random noise array

Return type:

jax.Array

step(key, field, force_updater=None)

Update the field by one time step.

Parameters:

• key (jax.random.PRNGKey) – Random number generator key

• field (Field) – The field to be updated

• force_updater (callable, optional) – A function that updates the force of fields

Returns:

The updated field

Return type:

Field

class openferro.integrator.llg.ConservativeLLSIBIntegrator(dt, gamma=None, max_iter=10, tol=1e-06)

Bases: Integrator

[Semi-implicit B (SIB) scheme. (J. Phys.: Condens. Matter 22 (2010) 176001) ] Adiabatic spin precession, i.e. Landau-Lifshitz equation of motion without dissipative damping term.

The equation of motion is

dM/dt = -gamma M x B

Let M[i] be the spin configuration at time step i, Y[i] be the auxiliary spin configuration at time step i+1, B(M) be the magnetic field of configuration M. The SIB scheme is

(step 1) Y[i] = M[i] - dt * gamma * (M[i]+Y[i])/2 x B(M[i])

(step 2) M[i+1] = M[i] - dt * gamma * (M[i]+M[i+1])/2 x B((M[i] + Y[i])/2)

Both equations are implicit, and can be solved iteratively through fixed-point iterations.

Parameters:

• dt (float) – Time step size

• gamma (float, optional) – Gyromagnetic ratio. If None, uses electron gyromagnetic ratio

• max_iter (int, optional) – Maximum number of iterations for fixed-point iterations

• tol (float, optional) – Tolerance for convergence. Convergence is declared if Average[|M_new - M_old|/Ms] < tol

_update_x(M, Ms, Y, B, step_size)

One iteration of (step 1) or (step 2) in the SIB scheme. Will be called iteratively through fixed-point iterations.

Parameters:

• M (jax.Array) – Current magnetization

• Ms (jax.Array) – Magnitude of magnetization

• Y (jax.Array) – Previous iteration result

• B (jax.Array) – Magnetic field

• step_size (float) – Time step size

Returns:

(Updated Y, normalized difference average)

Return type:

tuple

__init__(dt, gamma=None, max_iter=10, tol=1e-06)

step(field, force_updater)

Iteratively solve the implicit equations (step 1) and (step 2) in the SIB scheme.

Parameters:

• field (Field) – The field to be updated

• force_updater (callable) – A function that updates the force of fields

Returns:

The updated field

Return type:

Field

class openferro.integrator.llg.LLSIBIntegrator(dt, alpha, gamma=None, max_iter=10, tol=1e-06)

Bases: Integrator

Landau-Lifshitz equation of motion.

See Eriksson, Olle, et al. Atomistic spin dynamics: foundations and applications. Oxford university press, 2017, Sec.7.4.5 for details.

The equation of motion is

dM/dt = -gamma M x B - (gamma * alpha / |M|) * M x (M x B)

Here gamma=(gyromagnetic ratio)/ (1+alpha^2) is the renormalized gyromagnetic ratio for simulating LLG equation in Landau-Lifshitz form. Let M[i] be the spin configuration at time step i, Y[i] be the auxiliary spin configuration at time step i+1. Let B(M) be the effective magnetic field that includes also the dissipative damping term. The SIB scheme is

(step 1) Y[i] = M[i] - dt * gamma * (M[i]+Y[i])/2 x B(M[i])

(step 2) M[i+1] = M[i] - dt * gamma * (M[i]+M[i+1])/2 x B((M[i] + Y[i])/2)

Both equations are implicit, and can be solved iteratively through fixed-point iterations.

Parameters:

• dt (float) – Time step size

• alpha (float) – Gilbert damping constant

• gamma (float, optional) – Gyromagnetic ratio. If None, uses electron gyromagnetic ratio

• max_iter (int, optional) – Maximum number of iterations for fixed-point iterations

• tol (float, optional) – Tolerance for convergence. Convergence is declared if Average[|M_new - M_old|/Ms averaged over lattice] < tol

_update_x(M, Ms, Y, B, step_size)

One iteration of (step 1) or (step 2) in the SIB scheme. Will be called iteratively through fixed-point iterations.

Parameters:

• M (jax.Array) – Current magnetization

• Ms (jax.Array) – Magnitude of magnetization

• Y (jax.Array) – Previous iteration result

• B (jax.Array) – Magnetic field

• step_size (float) – Time step size

Returns:

(Updated Y, normalized difference average)

Return type:

tuple

_update_b(M, B, alpha, Ms)

Update the magnetic field B by adding the term ( alpha / |M|) * (M x B)

Parameters:

• M (jax.Array) – The magnetization (shape=(*, 3))

• B (jax.Array) – The magnetic field (shape=(*, 3))

• alpha (float) – The Gilbert damping constant

• Ms (jax.Array) – The magnitude of the magnetization (shape=(*,))

Returns:

The updated magnetic field (shape=(*, 3))

Return type:

jax.Array

__init__(dt, alpha, gamma=None, max_iter=10, tol=1e-06)

step(field, force_updater)

Iteratively solve the implicit equations (step 1) and (step 2) in the SIB scheme.

Parameters:

• field (Field) – The field to be updated

• force_updater (callable) – A function that updates the force of fields

Returns:

The updated field

Return type:

Field

class openferro.integrator.llg.LLSIBLangevinIntegrator(dt, temp, alpha, gamma=None, max_iter=10, tol=1e-06)

Bases: LLSIBIntegrator

Stochastic Landau-Lifshitz equation of motion.

See Eriksson, Olle, et al. Atomistic spin dynamics: foundations and applications. Oxford university press, 2017, Sec.7.4.5 for details.

The equation of motion is

dM/dt = -gamma M x (B + b) - (gamma * alpha / |M|) * M x (M x (B + b))

Here gamma=(gyromagnetic ratio)/ (1+alpha^2) is the renormalized gyromagnetic ratio for simulating LLG equation in Landau-Lifshitz form. b is the stochastic force <b_i_alpha(t) b_j_beta(s)> = 2 * D * delta_ij * delta_alpha_beta * delta(t-s) A steady Boltzmann state requires D= alpha/(1+alpha^2) * kbT / gamma / |m|

Let M[i] be the spin configuration at time step i, Y[i] be the auxiliary spin configuration at time step i+1. Let B(M[i]) be the effective magnetic field including also the dissipative damping term and the stochastic force. The SIB scheme is

(step 1) Y[i] = M[i] - dt * gamma * (M[i]+Y[i])/2 x B(M[i])

(step 2) M[i+1] = M[i] - dt * gamma * (M[i]+M[i+1])/2 x B((M[i] + Y[i])/2)

Both equations are implicit, and can be solved iteratively through fixed-point iterations.

__init__(dt, temp, alpha, gamma=None, max_iter=10, tol=1e-06)

get_noise(key, field)

step(key, field, force_updater=None)

Iteratively solve the implicit equations (step 1) and (step 2) in the SIB scheme.

Parameters:

• key (jax.random.PRNGKey) – Random number generator key

• field (Field) – The field to be updated

• force_updater (callable, optional) – A function that updates the force of fields

Returns:

The updated field with new spin configuration (shape=(*, 3))

Return type:

Field

Simulation

Classes which define the time evolution of physical systems.

Notes

This file is part of OpenFerro.

class openferro.simulation.Simulation(system)

Bases: object

The base class to define a simulation.

A simulation controls the time evolution of a system. This class does not implement any time integration algorithm. Each field object has its own integrator. The class only calls the step method of each integrator and controls output.

Parameters:

system (System) – The physical system to simulate

__init__(system)

clear_reporters()

add_thermo_reporter(file='thermo.log', log_interval=100, global_strain=False, excess_stress=False, volume=False, potential_energy=False, kinetic_energy=False, temperature=False)

Add a reporter to output global thermodynamic information.

Parameters:

• file (str, optional) – Output file name, by default ‘thermo.log’

• log_interval (int, optional) – Number of steps between outputs, by default 100

• global_strain (bool, optional) – Whether to output global strain, by default False

• excess_stress (bool, optional) – Whether to output excess stress, by default False

• volume (bool, optional) – Whether to output volume, by default False

• potential_energy (bool, optional) – Whether to output potential energy, by default False

• kinetic_energy (bool, optional) – Whether to output kinetic energy, by default False

• temperature (bool, optional) – Whether to output temperature, by default False

add_field_reporter(file_prefix, field_ID, log_interval=100, field_average=True, dump_field=False)

Add a reporter to dump the values of a given field.

Parameters:

• file_prefix (str) – Prefix for output files

• field_ID (str) – ID of field to report

• log_interval (int, optional) – Number of steps between outputs, by default 100

• field_average (bool, optional) – Whether to output field averages, by default True

• dump_field (bool, optional) – Whether to dump full field values, by default False

initialize_reporters()

Initialize all reporters.

remove_all_reporters()

reset_reporters()

Reset the counters of all reporters.

step_reporters()

Step all reporters.

init_velocity(mode='zero', temp=None)

_step()

Update the system by one time step. To be implemented by subclasses.

run()

Run the simulation for a given number of steps or until convergence. To be implemented by subclasses.

class openferro.simulation.MDMinimize(system, max_iter=100, tol=1e-05)

Bases: Simulation

Class for energy minimization using molecular dynamics.

Parameters:

• system (System) – The physical system to minimize

• max_iter (int, optional) – Maximum number of iterations, by default 100

• tol (float, optional) – Force tolerance for convergence, by default 1e-5

__init__(system, max_iter=100, tol=1e-05)

_step(variable_cell)

Update the field by one time step.

Parameters:

variable_cell (bool) – Whether to allow cell parameters to vary

run(variable_cell=True, pressure=None)

Run the minimization.

Parameters:

• variable_cell (bool, optional) – Whether to allow cell parameters to vary, by default True

• pressure (float, optional) – External pressure in bar, required if variable_cell=True

Raises:

ValueError – If pressure not specified for variable cell minimization If pressure specified for fixed cell minimization If integrator not set for any field

class openferro.simulation.SimulationNVE(system)

Bases: Simulation

Class for NVE (microcanonical) molecular dynamics simulation.

Parameters:

system (System) – The physical system to simulate

__init__(system)

_step(profile=False)

Update the field by one step.

Parameters:

profile (bool, optional) – Whether to profile timing, by default False

run(nsteps=1, profile=False)

Run the simulation.

Parameters:

• nsteps (int, optional) – Number of steps to run, by default 1

• profile (bool, optional) – Whether to profile timing, by default False

Raises:

ValueError – If integrator not set for any field

class openferro.simulation.SimulationNVTLangevin(system)

Bases: SimulationNVE

Class for NVT molecular dynamics simulation using Langevin dynamics.

Parameters:

system (System) – The physical system to simulate

__init__(system)

_step(keys, profile=False)

Update the field by one step.

Parameters:

• keys (array_like) – Random keys for Langevin dynamics

• profile (bool, optional) – Whether to profile timing, by default False

run(nsteps=1, profile=False)

Run the simulation.

Parameters:

• nsteps (int, optional) – Number of steps to run, by default 1

• profile (bool, optional) – Whether to profile timing, by default False

Raises:

ValueError – If integrator not set for any field

class openferro.simulation.SimulationNPTLangevin(system, pressure=0.0)

Bases: SimulationNVTLangevin

Class for NPT molecular dynamics simulation using Langevin dynamics.

Parameters:

• system (System) – The physical system to simulate

• pressure (float, optional) – External pressure in bar, by default 0.0

__init__(system, pressure=0.0)

Interaction (Abstract)

Classes which define the “interaction” between fields.

Each interaction is associated with a term in the Hamiltonian. Each interaction stores a function “self.energy_engine” that calculates the energy of the interaction and a function “force engine” that calculates the force of the interaction.

Only the energy engine is required. The force engine is optional. If the force engine is not set, the force will be calculated by automatic differentiation of the energy engine.

Notes

This file is part of OpenFerro.

class openferro.interaction.interaction_base(parameters=None)

Bases: object

The base class to specify the interaction between fields.

__init__(parameters=None)

set_parameters(parameters)

Set the parameters of the interaction.

Parameters:

parameters (array_like) – The parameters of the interaction

Raises:

ValueError – If parameters is not a numpy array, list, or jax array

get_parameters()

Get the parameters of the interaction.

Returns:

The parameters of the interaction

Return type:

jax.numpy.ndarray

set_energy_engine(energy_engine, enable_jit=True)

Set the energy engine of the interaction.

Parameters:

• energy_engine (callable) – The energy engine of the interaction. It should take the values of the fields as input and return the energy as output.

• enable_jit (bool, optional) – Whether to enable jit for the energy engine, by default True

calc_energy()

calc_force()

class openferro.interaction.self_interaction(field_ID, parameters=None)

Bases: interaction_base

A class to specify the self-interaction of a field.

Parameters:

• field_ID (str) – Identifier for the field

• parameters (array_like, optional) – Parameters for the interaction, by default None

__init__(field_ID, parameters=None)

create_force_engine(enable_jit=True)

Derive the force engine of the interaction from the energy engine through automatic differentiation.

Parameters:

enable_jit (bool, optional) – Whether to enable jit for the force engine, by default True

Raises:

ValueError – If energy engine is not set

calc_energy(field)

Calculate the energy of the interaction for a given field.

Parameters:

field (Field) – The field to calculate the energy

Returns:

The energy of the interaction

Return type:

float

calc_force(field)

Calculate the force of the interaction for a given field.

Parameters:

field (Field) – The field to calculate the force

Returns:

The gradient of the energy with respect to the field. It has the same shape as the field.

Return type:

jax.numpy.ndarray

class openferro.interaction.mutual_interaction(field_1_ID, field_2_ID, parameters=None)

Bases: interaction_base

A class to specify the mutual interaction between two fields.

Parameters:

• field_1_ID (str) – Identifier for the first field

• field_2_ID (str) – Identifier for the second field

• parameters (array_like, optional) – Parameters for the interaction, by default None

__init__(field_1_ID, field_2_ID, parameters=None)

create_force_engine(enable_jit=True)

Derive the force engine of the interaction from the energy engine through automatic differentiation.

Parameters:

enable_jit (bool, optional) – Whether to enable jit for the force engine, by default True

Raises:

ValueError – If energy engine is not set

calc_energy(field1, field2)

Calculate the energy of the interaction for a given pair of fields.

Parameters:

• field1 (Field) – The first field

• field2 (Field) – The second field

Returns:

The energy of the interaction

Return type:

float

calc_force(field1, field2)

Calculate the force of the interaction for a given pair of fields.

Parameters:

• field1 (Field) – The first field

• field2 (Field) – The second field

Returns:

The gradient of the energy with respect to the fields. It has the same shape as the fields.

Return type:

tuple of jax.numpy.ndarray

class openferro.interaction.triple_interaction(field_1_ID, field_2_ID, field_3_ID, parameters=None)

Bases: object

A class to specify the mutual interaction between three fields.

Parameters:

• field_1_ID (str) – Identifier for the first field

• field_2_ID (str) – Identifier for the second field

• field_3_ID (str) – Identifier for the third field

• parameters (array_like, optional) – Parameters for the interaction, by default None

__init__(field_1_ID, field_2_ID, field_3_ID, parameters=None)

create_force_engine(enable_jit=True)

Derive the force engine of the interaction from the energy engine through automatic differentiation.

Parameters:

enable_jit (bool, optional) – Whether to enable jit for the force engine, by default True

Raises:

ValueError – If energy engine is not set

calc_energy(field1, field2, field3)

Calculate the energy of the interaction for a given triple of fields.

Parameters:

• field1 (Field) – The first field

• field2 (Field) – The second field

• field3 (Field) – The third field

Returns:

The energy of the interaction

Return type:

float

calc_force(field1, field2, field3)

Calculate the force of the interaction for a given triple of fields.

Parameters:

• field1 (Field) – The first field

• field2 (Field) – The second field

• field3 (Field) – The third field

Returns:

The gradient of the energy with respect to the fields. It has the same shape as the fields.

Return type:

tuple of jax.numpy.ndarray

Utilities

Reporters

Classes for reporting simulation data.

This module contains classes for reporting thermodynamic quantities and field values during simulations.

class openferro.reporter.Thermo_Reporter(file='thermo.log', log_interval=100, global_strain=False, excess_stress=False, volume=True, potential_energy=False, kinetic_energy=False, temperature=False)

Bases: object

Reporter for thermodynamic quantities.

Parameters:

• file (str, optional) – Output file name, by default ‘thermo.log’

• log_interval (int, optional) – Number of steps between outputs, by default 100

• global_strain (bool, optional) – Whether to output global strain, by default False

• excess_stress (bool, optional) – Whether to output excess stress, by default False

• volume (bool, optional) – Whether to output volume, by default True

• potential_energy (bool, optional) – Whether to output potential energy, by default False

• kinetic_energy (bool, optional) – Whether to output kinetic energy, by default False

• temperature (bool, optional) – Whether to output temperature, by default False

__init__(file='thermo.log', log_interval=100, global_strain=False, excess_stress=False, volume=True, potential_energy=False, kinetic_energy=False, temperature=False)

initialize(system)

Initialize the reporter.

Parameters:

system (System) – The physical system to report on

step(system)

Record data for current step.

Parameters:

system (System) – The physical system to report on

class openferro.reporter.Field_Reporter(file_prefix, field_ID, log_interval=100, field_average=True, dump_field=False)

Bases: object

Reporter for field values.

Parameters:

• file_prefix (str) – Prefix for output files

• field_ID (str) – ID of field to report

• log_interval (int, optional) – Number of steps between outputs, by default 100

• field_average (bool, optional) – Whether to output field averages, by default True

• dump_field (bool, optional) – Whether to dump full field values, by default False

__init__(file_prefix, field_ID, log_interval=100, field_average=True, dump_field=False)

initialize(system)

Initialize the reporter.

Parameters:

system (System) – The physical system to report on

step(system)

Record data for current step.

Parameters:

system (System) – The physical system to report on

Units

Module containing classes that define the internal units.

Notes

This file is part of OpenFerro.

class openferro.units.Constants

Bases: object

Class containing fundamental physical constants and unit conversions. This class specifies the internal units of OpenFerro.

Angstrom = 1.0

length = 1.0

eV = 1.0

energy = 1.0

ps = 1.0

time = 1.0

elementary_charge = 1.0

charge = 1.0

Kelvin = 1.0

temperature = 1.0

derived units

kb = 8.6173303e-05

amu = 0.0001035

me = 5.677803e-08

epsilon0 = 0.005526349406

meter = 10000000000.0

cm = 100000000.0

mm = 10000000.0

um = 10000.0

nm = 10.0

pm = 0.001

bohr_radius = 0.52917721

second = 1000000000000.0

ms = 1000000000.0

us = 1000000.0

ns = 1000.0

fs = 0.001

bar = 6.24150912510413e-07

Coulomb = 6.241509074460763e+18

Ampere = 6241509.0744607635

Voltage = 1.0

V_Angstrom = 1.0

Joule = 6.241509e+18

Ry = 13.605693123

mRy = 0.013605693123000001

hbar = 0.0006582119569

hartree = 27.211386245988

we will use Bohr magneton based unit system for magnetic interactions. This is not entirely consistent with the other base units. For example, here mu0 will be in unit of eV * A^3 / muB^2, instead of eV * ps^2 / A * e^-2 Derivation: _speed_of_light = 2997924.58 # Speed of light in A/ps _mu0 = 1 / _speed_of_light**2 / epsilon0 # Vacuum permeability in eV*ps^2/A*e^-2 _muB = elementary_charge * hbar / 2 / me # Bohr magneton in e * A^2 / ps mu0 = _mu0 * _muB**2 / Angstrom**3 # Vacuum permeability in eV * A^3 / muB^2

muB = 1.0

mu0 = 0.0006764429231627333

Tesla = 5.788381806e-05

electron_g_factor = 2.00231930436

electron_gyro_ratio = 3042.058539608398

class openferro.units.AtomicUnits_to_InternalUnits

Bases: object

Class defining atomic units, given in terms of the internal units of OpenFerro.

time = 2.418884326505e-05

length = 0.52917721

energy = 27.211386245988

mass = 5.677803e-08

velocity = 21876.91259980993

force = 51.42206756407367

electric_dipole_moment = 0.52917721

Helper Functions

Utility functions.

openferro.utilities.SO3_rotation(B: Array, dt: float)

Calculate batched rotation matrix for proper rotation around magnetic field axis.

Parameters:

• B (jnp.ndarray) – The magnetic field with shape (*, 3)

• dt (float) – The time step

Returns:

The rotation matrix with shape (*, 3, 3)

Return type:

jnp.ndarray

Notes

Returns the batched rotation matrix (shape=(*, 3, 3)) of a proper rotation of a vector by an angle theta (shape=(*,)) around the axes u (shape=(*, 3)). See https://en.wikipedia.org/wiki/Rotation_matrix for details.

The angle and axes are specified by the magnetic field B and time step dt as: u = B / jnp.linalg.norm(B, axis=-1)[…, None] theta = jnp.linalg.norm(B, axis=-1) * dt

Parallelism

Classes for multi-GPU parallelism.

Notes

This file is part of OpenFerro.

class openferro.parallelism.DeviceMesh(devices=None, num_rows=None, num_cols=None)

Bases: object

__init__(devices=None, num_rows=None, num_cols=None)

Initialize the single-host multi-device parallelism. Get the mesh of the devices.

Parameters:

• devices (array-like, optional) – List of devices to use. If None, uses all available devices

• num_rows (int, optional) – Number of rows in device mesh. If None, automatically determined

• num_cols (int, optional) – Number of columns in device mesh. If None, automatically determined

Raises:

ValueError – If only one device is available If num_rows * num_cols does not match number of devices

partition_sharding()

Produce a NamedSharding object to distribute a value across devices, partitioning along the x and y axes.

Returns:

Sharding object for partitioning values across devices

Return type:

NamedSharding

replicate_sharding()

Produce a NamedSharding object to replicate a value across devices.

Used for broadcasting values that do not need to be partitioned.

Returns:

Sharding object for replicating values across devices

Return type:

NamedSharding