Ionic conductivity from ams.rkf trajectory

First, an ams.rkf trajectory with ions in the system needs to be calculated, as is done, for example, in the Ionic Conductivity Tutorial. Next, the ionic conductivity for the ions in the solution can be computed. The below script first identifies the ions in the molecular system, and guesses their charges (in the case of metals it uses a simple charge equilibration scheme). It then runs the AMS trajectory analysis tool to compute the ionic conductivity for these ions.

Example usage: (Download get_ionic_conductivity.py)

$AMSBIN/amspython get_ionic_conductivity.py /path/to/ams.rkf

To create an example ams.rkf file, a PLAMS script can be found here (Download NaClwater.py). This script is designed for efficiency. For a more reliable setup, follow the Ionic Conductivity Tutorial.

import sys
import os
from scm.plams import AMSJob, Molecule, Settings, Units, AMSAnalysisJob
from scm.plams.recipes.assign_charged_ions import assign_charged_ions


def main(filename):
    """
    The main body of the script
    """
    # Get the molecule
    job = AMSJob.load_external(filename)
    mol = job.molecule

    # Assign the charges to the ions
    mol = assign_charged_ions(mol, use_system_charges=False)

    # Read the temperature from KF
    job = AMSJob.load_external(filename)
    T = job.results.readrkf("MDResults", "MeanTemperature")
    kBT = Units.constants["Boltzmann"] * T
    print("Average temperaturs %f K" % (T))

    # Get the molecular system and extract the ions
    iontypes, ioncharges, nions, formulas = get_ions(mol)
    print("%8s %8s %10s" % ("Ion", "N", "Charge"))
    for k, indices in iontypes.items():
        print("%8s %8i %10.5f" % (formulas[k], len(indices), ioncharges[k]))

    # Compute diffusion coefficient for each ion
    diffusion_coeffs = {}
    for label, atoms in iontypes.items():
        s = Settings()
        s.input.Task = "MeanSquareDisplacement"
        s.input.TrajectoryInfo.Trajectory.KFFilename = filename
        atsettings = [iat + 1 for iat in atoms]
        s.input.MeanSquareDisplacement.Atoms.Atom = atsettings

        job = AMSAnalysisJob(settings=s)
        results = job.run()
        D = results._kf.read("Slope(1)", "Final")
        diffusion_coeffs[label] = D
        units = results._kf.read("Slope(1)", "Final(units)")
        print(formulas[label], D, units)

    # Compute the number density for each ion
    rho = {}
    for label, ni in nions.items():
        rho[label] = ni / mol.unit_cell_volume(unit="m")

    # Compute the ionic conductivity
    sigma = 0.0
    for label, D in diffusion_coeffs.items():
        q = ioncharges[label] * Units.constants["electron_charge"]
        s = q**2 * rho[label] * D / kBT
        sigma += s
    return sigma


def get_ions(mol):
    """
    Extract the ions and their charges from the atom properties
    """
    ions = {}
    for i, at in enumerate(mol.atoms):
        regions = []
        if "region" in at.properties:
            regions = at.properties.region
            if not isinstance(regions, set):
                regions = set([regions])
            regions = [s for s in regions if s[:3] == "ion"]
        if len(regions) == 0:
            continue
        name = regions[0]
        if not name in ions:
            ions[name] = []
        ions[name].append(i)

    # Extract the ion info
    iontypes = {}
    ioncharges = {}
    nions = {}
    formulas = {}
    for name, atoms in ions.items():
        charge = 0
        for iat in atoms:
            if "analysis" in mol.atoms[iat].properties:
                if "charge" in mol.atoms[iat].properties.analysis:
                    charge += mol.atoms[iat].properties.analysis.charge
                else:
                    raise PlamsError("Not all charges present in system file")
            elif "forcefield" in mol.atoms[iat].properties:
                if "charge" in mol.atoms[iat].properties.forcefield:
                    charge += mol.atoms[iat].properties.forcefield.charge
                else:
                    raise PlamsError("Not all charges present in system file")
        typename = name.split("_")[0]
        if not typename in iontypes:
            ion = mol.get_fragment(atoms)
            formulas[typename] = ion.get_formula()
            ioncharges[typename] = charge
            iontypes[typename] = []
            nions[typename] = 0
        iontypes[typename] += atoms
        nions[typename] += 1

    return iontypes, ioncharges, nions, formulas


if __name__ == "__main__":
    if len(sys.argv) < 2:
        print("Usage: amspython get_ionic_conductivity.py path/to/ams.rkf")
        sys.exit(0)
    filename = os.path.abspath(sys.argv[1])
    sigma = main(filename)
    print("Ionic conductivity: %20.10e Siemens/m" % (sigma))