Python Scripts

This page describes the collection of Python sample scripts contained in the Python sample scripts folder.

Run a VASP relaxation using the internal quasi-Newton optimization algorithm (RMM-DIIS) in VASP

samples/python/run_vasp.py

from pathlib import Path

from ase import io
from ase.calculators.vasp import Vasp
from numpy.linalg import norm

# Replace in.traj with the name of your structure file
atoms = io.read("in.traj")

# see https://www.vasp.at/wiki/index.php/Category:INCAR_tag
# for details on what each of these keywords mean
# if not found, check https://wiki.fysik.dtu.dk/ase/ase/calculators/vasp.html#module-ase.calculators.vasp
calc = Vasp(
    algo="Fast",
    encut=450,
    gga="PE",
    gamma=False,
    ibrion=1,
    isif=2,
    ismear=0,
    ncore=4,
    nelm=60,
    nsw=100,
    prec="Accurate",
    sigma=0.1,
    kpts=(4, 4, 1),
)

atoms.calc = calc

e = atoms.get_potential_energy()
f = norm(max(atoms.get_forces(), key=norm))

# Write final structure to file
atoms.write("final.traj")

# Print final energy and max force to standard output
print(f"final energy {e}")
print(f"max force: {f}")

# Write final energy to file
with Path("final.e").open(mode="x", encoding="utf-8") as file:
    file.write(f"{e}\n")

Reminder

Replace "in.traj" with the name of your structure file.

Run a mixed-basis Gaussian relaxation using the internal Berny optimization algorithm in Gaussian

samples/python/run_gaussian.py

from pathlib import Path

from ase import io
from ase.calculators.gaussian import Gaussian
from numpy.linalg import norm

# Replace in.traj with the name of your structure file
atoms = io.read("in.traj")

# see https://gaussian.com/keywords/
# for details on what each of these keywords mean
# if not found, check https://wiki.fysik.dtu.dk/ase/ase/calculators/gaussian.html
calc = Gaussian(
    command="g16 < Gaussian.com > Gaussian.log",
    label="Gaussian",
    chk="Gaussian",
    save=None,
    xc="pbe0",
    basis="gen",
    basis_set="-H -N -C -O\nDef2SVP\n****\n-Cu -Fe\nDef2TZVP\n****\n",
    charge=0,
    mult=3,
    scf="QC,MaxCycle=100,IntRep",
    opt="Tight,MaxCycles=250,CalcFC",
    mem="40GB",
)

atoms.calc = calc

e = atoms.get_potential_energy()
f = norm(max(atoms.get_forces(), key=norm))

# Write final structure to file
atoms.write("final.traj")

# Print final energy and max force to standard output
print(f"final energy {e}")
print(f"max force: {f}")

# Write final energy to file
with Path("final.e").open(mode="x", encoding="utf-8") as file:
    file.write(f"{e}\n")

Reminder

Replace "in.traj" with the name of your structure file, and read this article about selecting which executable to pass to the command keyword argument to the Gaussian constructor. Note that the syntax is slightly different (g16 < Gaussian.com vs. G16 Gaussian.com).

Run a VASP relaxation using the BFGSLineSearch optimization routine in ASE

samples/python/run_ase.py

from pathlib import Path

from ase import io
from ase.calculators.vasp import Vasp
from ase.optimize.bfgslinesearch import BFGSLineSearch
from numpy.linalg import norm

# Replace in.traj with the name of your structure file
atoms = io.read("in.traj")

# see https://www.vasp.at/wiki/index.php/Category:INCAR_tag
# for details on what each of these keywords mean
# if not found, check https://wiki.fysik.dtu.dk/ase/ase/calculators/vasp.html#module-ase.calculators.vasp
calc = Vasp(
    algo="Fast",
    gga="PE",
    gamma=False,
    ibrion=1,
    isif=2,
    ismear=0,
    nelm=60,
    # Note that nsw must equal zero in order to properly  use the ASE optimizer
    nsw=0,
    prec="Accurate",
    kpts=(4, 4, 1),
)

atoms.calc = calc

# Create the ASE optimizer and run the optimization
dyn = BFGSLineSearch(
    atoms,
    logfile="BFGS_output.log",
    trajectory="relax.traj",
    restart="relax.pckl",
)
dyn.run(fmax=0.01, steps=150)

e = atoms.get_potential_energy()
f = norm(max(atoms.get_forces(), key=norm))

# Write final structure to file
atoms.write("final.traj")

# Print final energy and max force to standard output
print(f"final energy {e}")
print(f"max force: {f}")

# Write final energy to file
with Path("final.e").open(mode="x", encoding="utf-8") as file:
    file.write(f"{e}\n")

Reminder

Replace "in.traj" with the name of your structure file.

Perform a vibrational calculation using VASP

samples/python/run_vib.py

from pathlib import Path

from ase import io
from ase.calculators.vasp import Vasp
from ase.constraints import FixAtoms
from ase.vibrations.vibrations import Vibrations
from numpy.linalg import norm

# Replace in.traj with the name of your structure file
atoms = io.read("in.traj")

# see https://www.vasp.at/wiki/index.php/Category:INCAR_tag
# for details on what each of these keywords mean
# if not found, check https://wiki.fysik.dtu.dk/ase/ase/calculators/vasp.html#module-ase.calculators.vasp
calc = Vasp(
    algo="Fast",
    encut=450,
    gga="PE",
    gamma=False,
    ibrion=1,
    isif=2,
    ismear=0,
    ncore=4,
    nelm=60,
    nsw=100,
    prec="Accurate",
    sigma=0.1,
    kpts=(4, 4, 1),
)

atoms.calc = calc

e = atoms.get_potential_energy()
f = norm(max(atoms.get_forces(), key=norm))

# Write final structure to file
atoms.write("final.traj")

# Print final energy and max force to standard output
print(f"final energy {e}")
print(f"max force: {f}")

# Write final energy to file
with Path("final.e").open(mode="x", encoding="utf-8") as file:
    file.write(f"{e}\n")

print("Running vibration calculation")
forces = atoms.get_forces(apply_constraint=True)
atoms.write("final.traj")

if forces is None:
    msg = "Unable to perform vibrational calculation due to missing forces"
    raise RuntimeError(msg)

print(">>> BEGIN print full force")

for force, atom in zip(forces, atoms, strict=True):
    print("%s %s: %s, %s, %s", atom.index, atom.symbol, *force)

print("norm of force: %s", norm(forces))
print("max force: %s", norm(max(forces, key=norm)))
print("<<< END print full force")

fixed = []
if atoms.constraints:
    for constraint in atoms.constraints:
        if isinstance(constraint, FixAtoms):
            fixed = constraint.index
            break

print("Fixed indices: %s", fixed)

to_vibrate = [atom.index for atom in atoms if atom.index not in fixed]
print("Indices to vibrate: %s", to_vibrate)

vib = Vibrations(
    atoms=atoms, nfree=4, name="vib", delta=0.015, indices=to_vibrate
)
vib.run()

with Path("vib.txt").open(mode="w", encoding="utf-8") as file:
    vib.summary(log=file)

for mode in range(len(to_vibrate) * 3):
    vib.write_mode(mode)

zpe = vib.get_zero_point_energy()

print("Successfully ran vibration calculation. ZPE: %s", zpe)

Run a VASP relaxation using ccu

samples/python/run_ccu.py

import logging

from ase.calculators.vasp import Vasp
import ase.io
from ccu.relaxation import run_relaxation

logging.basicConfig(level=logging.DEBUG)

atoms = ase.io.read("in.traj")

calc = Vasp(
    gga="PE",
    gamma=False,
    ibrion=1,
    isif=2,
    ismear=0,
    kpts=(1, 1, 1),
    nelm=60,
    nsw=100,
    prec="Accurate",
)

atoms.calc = calc
run_relaxation(atoms, run_bader=False, run_chargemol=False)

Perform a vibrational calculation using VASP and ccu

samples/python/run_vib_ccu.py

import logging

from ase.calculators.vasp import Vasp
import ase.io
from ccu.thermo.vibration import run_vibration

logging.basicConfig(level=logging.DEBUG)

# Replace in.traj with the name of your structure file
atoms = ase.io.read("in.traj")

# see https://www.vasp.at/wiki/index.php/Category:INCAR_tag
# for details on what each of these keywords mean
# if not found, check https://wiki.fysik.dtu.dk/ase/ase/calculators/vasp.html#module-ase.calculators.vasp
calc = Vasp(
    algo="Fast",
    encut=450,
    gga="PE",
    gamma=False,
    ibrion=1,
    isif=2,
    ismear=0,
    ncore=4,
    nelm=60,
    nsw=100,
    prec="Accurate",
    sigma=0.1,
    kpts=(4, 4, 1),
)

atoms.calc = calc
run_vibration(atoms)

Perform an infrared calculation using VASP and ccu

samples/python/run_infrared_ccu.py

import logging

from ase.calculators.vasp import Vasp
import ase.io
from ccu.thermo.infrared import run_infrared

logging.basicConfig(level=logging.DEBUG)

# Replace in.traj with the name of your structure file
atoms = ase.io.read("in.traj")

# see https://www.vasp.at/wiki/index.php/Category:INCAR_tag
# for details on what each of these keywords mean
# if not found, check https://wiki.fysik.dtu.dk/ase/ase/calculators/vasp.html#module-ase.calculators.vasp
calc = Vasp(
    algo="Fast",
    encut=450,
    gga="PE",
    gamma=False,
    ibrion=1,
    isif=2,
    ismear=0,
    ncore=4,
    nelm=60,
    nsw=100,
    prec="Accurate",
    sigma=0.1,
    kpts=(4, 4, 1),
)

atoms.calc = calc
run_infrared(atoms)

ccu is a set of tools for computational chemistry workflows. In particular, run_relaxation [run_vibration] and [run_infrared] function wrap the ase.atoms.Atoms.get_potential_energy(), ase.vibrations.vibrations.Vibration, and se.vibrations.infrared.Infrared functions, respectively and handle the logging and archiving of a calculation's final results.

Reminder

Replace "in.traj" with the name of your structure file.

Run an AIMD calculation with Quantum Espresso

samples/python/run_espresso.py

from pathlib import Path

from ase import io
from ase.calculators.espresso import Espresso
from ase.calculators.espresso import EspressoProfile
from numpy.linalg import norm

# Replace in.traj with the name of your structure file
atoms = io.read("in.traj")

# see https://www.quantum-espresso.org/Doc/INPUT_PW.html
# for details on what each of these keywords mean
# if not found, check https://wiki.fysik.dtu.dk/ase/ase/calculators/espresso.html
input_data = {
    "calculation": "md",
    "verbosity": "low",
    "restart_mode": "from_scratch",
    "tstress": False,
    "tprnfor": False,
    # This should be set to less than the job time limit
    "max_seconds": 75000,
    "occupations": "smearing",
    "smearing": "gaussian",
    "degauss": 0.1,
    "vdw_corr": "grimme-d3",
    "dftd3_version": 4,
}

# Ensure that you define an entry for each atom type in your structure
pseudopotentials = {
    "C": "C.pbe-n-kjpaw_psl.1.0.0.UPF",
    "Co": "Co_pbe_v1.2.uspp.F.UPF",
    "S": "s_pbe_v1.4.uspp.F.UPF",
}

profile = EspressoProfile(
    # Replace the path in pseudo_dir with the path to your QE pseudopotential
    # folder
    pseudo_dir="/home/ugognw/projects/def-samiras/ugognw/software_support/espresso/SSSP_1.3.0_PBE_precision",
    command="mpiexec pw.x -nband 2 -ntg 2",
)

calc = Espresso(
    input_data=input_data,
    pseudopotentials=pseudopotentials,
    profile=profile,
)

atoms.calc = calc

e = atoms.get_potential_energy()
f = norm(max(atoms.get_forces(), key=norm))

# Write final structure to file
atoms.write("final.traj")

# Print final energy and max force to standard output
print(f"final energy {e}")
print(f"max force: {f}")

# Write final energy to file
with Path("final.e").open(mode="x", encoding="utf-8") as file:
    file.write(f"{e}\n")

Reminder

Replace "in.traj" with the name of your structure file, and don't forget to change the string passed to the pseudo_dir keyword argument of the EspressoProfile constructor to the absolute path to where your Quantum Espresso pseudopotentials are located.