Volume Scan with Birch-Murnaghan Equation of State¶

Run a lattice PES scan for diamond with ReaxFF, fit the energy-volume curve to a Birch-Murnaghan equation of state, and rescale the cell to the fitted equilibrium volume.

Downloads: Notebook | Script ?

Requires: AMS2027 or later

../_images/volume_scan_birch_murnaghan_7_0_548729e6.png

Initial structure¶

To get the initial diamond structure in AMSinput:

Builders -> Crystal Structure -> Cubic -> Diamond
Click “Ok”
File -> “Export Coordinates (System)” -> “.in…”
Copy the contents of the resulting .in file

from scm.base import ChemicalSystem
from scm.plams import view

diamond = ChemicalSystem("""
System
    Atoms
        C -0.44625 -0.44625 -0.44625
        C 0.44625 0.44625 0.44625
    End
    Lattice
        0.0 1.785 1.785
        1.785 0.0 1.785
        1.785 1.785 0.0
    End
End
""")

view(diamond, direction="tilt_x", guess_bonds=True)

Run a volume scan for the diamond structure¶

from scm.inputs import AMS, ReaxFF
from scm.plams import AMSJob

sett = AMS(Task="PESScan")
sett.PESScan.ScanCoordinate = [{"CellVolumeScalingRange": [0.8, 1.3], "nPoints": 11}]
sett.Engine = ReaxFF(ForceField="C.ff")

job = AMSJob(name="DiamondVolumeScan", settings=sett, molecule=diamond)
job.run();

[08.10|14:50:31] JOB DiamondVolumeScan STARTED
[08.10|14:50:31] JOB DiamondVolumeScan RUNNING
[08.10|14:50:32] JOB DiamondVolumeScan FINISHED
[08.10|14:50:32] JOB DiamondVolumeScan SUCCESSFUL

Extract the results¶

Obtain the volumes in Å3
Obtain the energies in eV
Plot the results and fit with the Birch-Murnaghan EoS

from scm.utils.conversions import plams_molecule_to_chemsys
from scm.base import Units

from ase.eos import EquationOfState
import matplotlib.pyplot as plt

res = job.results.get_pesscan_results(molecules=True)

# Obtain volumes
volumes = [plams_molecule_to_chemsys(x).lattice.get_volume() for x in res["Molecules"]]

# Obtain energies
energy_unit_conversion = Units.conversion_factor("Ha", "eV")
energies = [energy * energy_unit_conversion for energy in res["PES"]]

# Fit and plot
eos = EquationOfState(volumes, energies, eos="birchmurnaghan")
v0, e0, B = eos.fit()
eos.plot(show=True);

print(f"{e0=:8.3f} eV\n{v0=:8.3f} Å\u00b3\n{ B=:8.3f} GPa")

e0= -15.150 eV
v0=  11.463 Å³
 B=   4.439 GPa

Scale cell to match minimum volume¶

set_density applies uniform scaling.

diamond_scaled = diamond.copy()
diamond_scaled.set_density(diamond.total_mass() / v0)

print(diamond_scaled)

System
   Atoms
      C -0.44739431927580825 -0.44739431927580825 -0.44739431927580825
      C 0.44739431927580825 0.44739431927580825 0.44739431927580825
   End
... output trimmed ....
      0 1.7895772771032332 1.7895772771032332
      1.7895772771032332 0 1.7895772771032332
      1.7895772771032332 1.7895772771032332 0
   End
End

print(f"{diamond.lattice.get_volume()=:.3f}")
print(f"{diamond_scaled.lattice.get_volume()=:.3f}")

diamond.lattice.get_volume()=11.375
diamond_scaled.lattice.get_volume()=11.463

Python Script¶

#!/usr/bin/env python
# coding: utf-8

# ## Initial structure

# To get the initial diamond structure in AMSinput:
# 
# 1. Builders -> Crystal Structure -> Cubic -> Diamond
# 2. Click "Ok"
# 3. File -> "Export Coordinates (System)" -> ".in..."
# 4. Copy the contents  of the resulting .in file

from scm.base import ChemicalSystem
from scm.plams import view

diamond = ChemicalSystem("""
System
    Atoms
        C -0.44625 -0.44625 -0.44625 
        C 0.44625 0.44625 0.44625 
    End
    Lattice
        0.0 1.785 1.785
        1.785 0.0 1.785
        1.785 1.785 0.0
    End
End
""")


view(diamond, direction="tilt_x", guess_bonds=True, picture_path="picture1.png")


# ## Run a volume scan for the diamond structure

from scm.inputs import AMS, ReaxFF
from scm.plams import AMSJob

sett = AMS(Task="PESScan")
sett.PESScan.ScanCoordinate = [{"CellVolumeScalingRange": [0.8, 1.3], "nPoints": 11}]
sett.Engine = ReaxFF(ForceField="C.ff")

job = AMSJob(name="DiamondVolumeScan", settings=sett, molecule=diamond)
job.run();


# ## Extract the results
# 
#  - Obtain the volumes in Å<sup>3</sup>
#  - Obtain the energies in eV
#  - Plot the results and fit with the Birch-Murnaghan EoS

from scm.utils.conversions import plams_molecule_to_chemsys
from scm.base import Units

from ase.eos import EquationOfState
import matplotlib.pyplot as plt

res = job.results.get_pesscan_results(molecules=True)

# Obtain volumes
volumes = [plams_molecule_to_chemsys(x).lattice.get_volume() for x in res["Molecules"]]

# Obtain energies
energy_unit_conversion = Units.conversion_factor("Ha", "eV")
energies = [energy * energy_unit_conversion for energy in res["PES"]]

# Fit and plot
eos = EquationOfState(volumes, energies, eos="birchmurnaghan")
v0, e0, B = eos.fit()
eos.plot(show=True);


print(f"{e0=:8.3f} eV\n{v0=:8.3f} Å\u00b3\n{ B=:8.3f} GPa")


# ## Scale cell to match minimum volume
# 
# `set_density` applies uniform scaling.

diamond_scaled = diamond.copy()
diamond_scaled.set_density(diamond.total_mass() / v0)


print(diamond_scaled)


print(f"{diamond.lattice.get_volume()=:.3f}")
print(f"{diamond_scaled.lattice.get_volume()=:.3f}")

Volume Scan with Birch-Murnaghan Equation of State¶

Initial structure¶

Run a volume scan for the diamond structure¶

Extract the results¶

Scale cell to match minimum volume¶

See also¶

Related examples¶

Related documentation¶

Python Script¶