Calculating Elastic Constants

When you pull on a rubber band and then let go it does what? Snaps back to its original shape or length! This is elasticity, an applied displacement/strain results in a stored energy or restoring force/stress which upon release brings the material back to its original state (mostly¹). This is the theory of linear elasticity, meaning that the resulting stress is linear in response to the applied strain along with some coefficient which is material dependent, more specifically:

\begin{equation}\label{eq:linear_elasticity} \sigma_{ij} = C_{ijkl} \epsilon_{kl} \end{equation}

with the indices corresponding to the tensor Einstein notation and $\sigma$, $\epsilon$, and $C$ begin the stress, strain, and stiffness tensors²³. Its pretty straight forward as eq. \ref{eq:linear_elasticity} is linear as mentioned. There is of course non-linear elasticity which maybe important when dealing with materials that exhibit hyper-elasticity or coupling physics (e.g., multi-ferroics).

The measurement of elastic constants or the stiffness tensor components is a farily involved process since it requires preparation of single crystals and subsequent application of a force or strain along specific axes and facets of the crystal. However, because of conservation principles, the number of elastic constants is constrained by the crystal space group symmetry, and therefore for many crystals (e.g., Cubic crystals) the number of measurements is significantly reduced.

Now a days, computational techniques enable determination of all the elements in the elastic constant tensor. Since most of the focus is usually on linear elasticity, its very simple to perform these types of calculations using atomic displacements and then storing the resulting energy and forces. From this data, one can fit a the linear response and find the slope corresponds to the specific element of the elastic constant tensor, as shown below:

Example of MgO from the Elastic ASE extension package. Shows the linear response of stress vs. strain and from the fit one gets the elements of the elastic constant tensor. Note that they don't pass through the origin, which is probably because the residual stress/pressure of the relaxed VASP cell hasn't been accounted for.

If we have a cubic crystal, then by symmetry there are only 3 elastic constants in the stiffness tensor that we have to measure or calculate:

\begin{equation} \label{eq:elastic_tensor} \begin{pmatrix} C_{11} & C_{12} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{12} & C_{11} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{44} \end{pmatrix} \end{equation}.

There are some numerical details that makes determining elements of $C_{ijkl}$ a bit more simple to do. In particular, eq. \ref{eq:linear_elasticity} is transformed into a matrix representation where the strain matrix elements (i.e., strain vectors) encapsulate the respective symmetry, as a result one gets a $n \times 6$ set of linear equations for stress and strain because of different strain vectors data points $n$ one might perform. The shape of the strain matrix will depend on the number of independent elastic constants (i.e., symmetry), thus for a cubic crystal it would be $6 \times 3$ matrix multiplied by a vector of length $3$ to give the resulting symmetrized stress in vector form. However, we are trying to find the unique (i.e., generalized eigen-decomposition) vector corresponding to the elastic constants, thus we need to use a numerical technique that iterates over the different data points for the strain and stress matrices we have of a given crystal symmetry. To do this one can solve least-squares using SVD. More details on implementation and theory can be found in refs. [1-2].

Using ASE Elastic module

Much of the heavy lifting for these types of calculations is in preparing the crystal symmetry details and prescribing the deformations. This is where the usefulness of the Elastic ASE module [1] comes in handy. With as little as three lines, you can perform the calculations to get the elastic constants:

# Import ASE, create structure, set calculator


from elastic import get_elementary_deformations
from elastic.elastic import get_cij_order
from elastic import get_elastic_tensor


systems = get_elementary_deformations(structure, n=10, d=0.1)
result = SerialCalculate(systems, calculator)
ordering = get_cij_order(structure)


cij, _ = get_elastic_tensor(structure, systems=result)
cij /= units.GPa
elastic_constants = dict(zip(ordering, cij))

The structure and calculator variables correspond to ASE objects, so follow your respective steps.

Note

The Elastic ASE module was developed predominetly with VASP, Aims, and Siesta electronic structure codes in mind. However, I don't see anything that precludes using other calculators and have written a class for the LAMMPSlib that seems to work without any issues, other than it has to be run in serial mode (i.e. runs a single displacement calc at a time, but is fast for classical potentials). I opened an issue on GitHub that shows some example code.

The systems variable stores the deformations for the respective crystal symmetry. The variable n indicate how many instances , for an elementary deformation, to generate along ranges of -d to d which are the min and max deformation percentages about various axes and angles. The get_cij_order just provides the string value list of elastic components for a given symmetry, so that you can combine that with the result from the get_elastic_tensor function. The last not to keep in mind is that the units of cij will depend on your calculator setting in ASE.

With these simple lines of python code you can calculate the elastic constants using strain-stress linear response. Numerically, it is fairly accurate although there are other approaches based on strain-energy formalism that can be used.

Footnotes

---

1. In usual application of elasticity theory one assumes strict linear response and small displacements and thus the material returns to the equilibrium state upon release. Anything that deviates is treated in a different manner, for example, viscoelasticity or plasticity. ↩

2. Tensors are a convenient mathematical way to represent relationships/transformations of multi-dimensional quantities. For example its can become cumbersome to try and write what the force vector along a face of a Triclinic crystal would be given some displacement along a direction. Tensors allow for the "contraction" along shared indices. ↩

3. I'm going to use stiffness tensor and elastic constant tensor interchangeably. The inverse is to refer to the compliance tensor, which is the response to a stress rather than strain, but these are just inverse relations. ↩

References

[1] Paweł T. Jochym, Codacy Badger, Elastic: ASE Module, (2018). https://doi.org/10.5281/ZENODO.1254570. [2] L.D. Landau, E. Lifshitz, Theory of elasticity, 3. Engl. ed., rev.enlarged, Elsevier, Butterworth-Heinemann, Amsterdam Heidelberg, 2009. URL.

Reuse and Attribution

Custodians of AI: Ferrari Turbocharger or Brembo Brakes?

Right before Thanksgiving, the AI world seemed to unravel with the bizarre events at OpenAI: unexpected removal of the CEO and then his abrupt return. From my perspective, this drama sparked a monumental debate among the online AI community on platforms like LinkedIn, Twitter, and Reddit, framing it as a colossal clash of titans over the progression towards artificial general intelligence. I want to provide a quick overview of my understanding into this intriguing subculture that's emerged around this controversy.

The Rise of the Accelerationist

The term "effective accelerationism", abbreviated online as e/acc, seems to have gained a lot of prominence in the last 6 months with the discussion around LLMs and artificial general intelligence. Initially, I correlated it with "effective altruism" – a concept centered on maximizing societal good¹. However, effective accelerationism diverges from altruism in very specific way by proposing rapid technological advancements as the vehicle for societal betterment [1,2]². This approach entails moving at swift, often unchecked, velocities towards technological advancement. In my view, this is the embodying essence of engineering under the guise of a grander philosophical vision for a Utopian society. Essentially, e/acc is nothing more than engineering infused with some deeper philosophical aspirations. At least this is my take on it.

The Decelerationist Stance

Contrasting the accelerationists are the decelerationists, or decels³, who seem to be primarily AI researchers concerned with AI system safety [3]. Their apprehension seems two-fold, but it is unclear to me which is the most concerning. First, the fear of autonomous AI systems gaining unpredictable agency and, secondly, the potential misuse of AI by humans, like using AI-generated text or images to spread damaging misinformation. The decelerationists advocate for a thorough, safety-first approach above all. Pretty sure this is the mindset of most Ph.D. scientist – ever-questioning, often sidetracked by hypotheticals, and meticulously slow in deriving scientific conclusions. This has great value though! Without this mindset you wouldn't have the computer your reading this on.

However, me casting this as a binary conflict is far too simplistic. The quest for technological advancement isn't a zero-sum game, if the technologist prevail we don't for certainty if AGI or super-intelligence that eliminates us. If AI safety experts, lobby for regulation that slows the pace of innovation, humanity doesn't become an island of despair. Effective accelerationists just engineers and technologists who are visionaries eager to propel humanity into a future reminiscent of the Sci-Fi novels from their childhood. They probably pursue this with a sense of responsibility, not some mad scientist bent on inventing life, albeit tinged with excitement. Conversely, decelerationists, while cautious, are equally committed to leveraging technology for a brighter future. Their apprehensions often stem from observed disappointments from observing how platforms like TikTok can manipulate human behavior and interactions through algorithmic social engineering.

My "Genius" Thoughts, 🤪

Whats my takeaway from the e/acc and decel subculture? It's more than a simplistic 'us versus them' narrative. The online activities that label it as such is reductive and unproductive. It's a discourse rooted in good intentions, where decelerists aren't merely obstructionists, and accelerationists aren't recklessly charging forward without consideration. At its core, advances in technology are just about responsible engineering – a balanced approach to innovation and caution.

Footnotes

---

1. I have a pretty strong bias against the effective alturism movement. It might be a "good" guiding principal, but in actuation, it will never work. ↩

2. Although I used this reference I strongly diasgree with the practice of "doxxing" simply because the article author/editors viewed it as a "public good". Its interesting because I took a bit of personal offense because I'm very familiar with the quantum computing work by the person who was doxxed; he has some really nice quantum machine learning papers. But of course thats my problem. ↩

3. I'm curious if this is a slight at this group of people, it seems, uncannily close to incels ↩

References

[1] E. Baker-White, Who Is @BasedBeffJezos, The Leader Of The Tech Elite’s ‘E/Acc’ Movement?, Forbes. Article URL (accessed December 5, 2023).

[2] Effective accelerationism, Wikipedia. (2023). Wikipedia URL (accessed December 5, 2023).

[3] Pause Giant AI Experiments: An Open Letter, Future of Life Institute. https://futureoflife.org/open-letter/pause-giant-ai-experiments/ (accessed December 5, 2023).

Reuse and Attribution

Unofficial way to use Phonopy with ASE

If you've ever done any phonon calculations with VASP, then you are most certainly familiar with phonopy. In many scenarios you may want to use interatomic potentials for the forces instead of ab-initio¹ calculated forces. This allows you to more quickly converge and calculate ground-state bandstructures and thermal quantities. You can even look at finite temperature effects with lammps using fix phonon. Phonopy does have a good amount of support for several DFT² and LAMMPS.

I recently have been using atomic simulaton environment (ASE) to leverage some of its calculator interfaces. While ASE does have built in finite-difference method support for phonon calculations, phonopy is more comprehensive. So I wanted to explore how much effort it would take to get the two to "talk", turns out its extremly straightfoward.

The Workflow

First things first, you'll need to run pip install ase and pip install phonopy, its also ideal to install the Brillouin zone path library via pip install seekpath. Now, let me create a structure with ASE and define a calculator, for simplicity I'm going to use the effective medium potential (EMT) calculator thats built into ASE.

from ase.build import bulk
from ase.calculators.emt import EMT
calculator = EMT()
structure = bulk('Cu', 'fcc', a=3.6)

Now to import the phonopy python API tools:

from phonopy import Phonopy
from phonopy.interface.calculator import get_force_sets
from phonopy.structure.atoms import PhonopyAtoms

The import thing to notice is we need to use PhonopyAtoms to create the supercells and displacements. At first this may seem like a big hurdle, in that you'll need a lot of wrapper code, but this is not the case. Simply do:

phnpy_struct = PhonopyAtoms(
    symbols=structure.get_chemical_symbols(),
    positions=structure.get_positions(),
    cell=structure.get_cell(),
    )

Now we construct our phonopy object:

phonons = Phonopy(
    phnpy_struct,
    supercell_matrix=[[5, 0, 0],
                      [0, 5, 0],
                      [0, 0, 5]],
    )

Then we need to specify the displacement magnitude that is used when using the finite-difference method.

phonon.generate_displacements(distance=0.5)

Note

One thing I've noticed is that the distance values to converge the phonon calculations are significantly larger than I've used in the past. Usually displacements of 0.01-0.05 angstroms, but I'm noticing 0.5-0.8 angstroms are need to converge! Seems concerning since we need to remain within the harmonic regime.

Okay straightforward, now what about the forces? How do we get them into our phonopy object phonons. Turns out this is also relatively simple:

import numpy as np
sets_of_forces = []
supercells = phonons.get_supercells_with_displacements()
for i,d in enumerate(supercells):
    # Convert back to ASE
    d_ase = Atoms(symbols=d.get_chemical_symbols(),
        positions=d.get_positions(),
        cell=d.get_cell()
        )
    d_ase.set_calculator(calculator)
    forces = d_ase.get_forces()
    sets_of_forces.append(forces)

Now that we have the forces, fortunately, all we need to do to get the force constants to calculate the phonon properties is:

sets_of_forces = np.array(sets_of_forces)
phonons.forces = sets_of_forces
phonons.produce_force_constants()
phonons.save()

The last command saves the details for phonopy which I believe you can use the outputs to run the standalone phonopy command-line tool. Now we can view the phonon bandstructure:

phonons.phonons.auto_band_structure(plot=True)

which gives:

Phonon bandstructure of Cu FCC using ASE and Phonopy together.

Using the phonopy python api you should also be able to calculate the DOS and thermal properties.

Footnotes

---

1. If the forces on atoms are calculated using the Hellmann-Feynman theorem, then they are referred to as first-principles or ab-initio because they arise from the expectation value of the Hamiltonian about some continous parameter, ex., $\mathbf{F}_i = -\frac{dE_i}{d\mathbf{r}_i} = \langle \psi_{\mathbf{r_i}} \lvert \frac{d\hat{H}_{\mathbf{r}_i}}{{d\mathbf{r}_i}} \rvert \psi_{\mathbf{r_i}} \rangle$. ↩

2. Density functional theory is the primary compute for most high-quality phonon calculations. ↩

Reuse and Attribution