Skip to content

Convert from Amber restart to Geometry.dat file for AIMS🔗

scripts.convert_to_geomdat 🔗

Written by Alice Walker and Ruibin Liang

This script converts from Amber restart format to AIMS Geometry.dat files. It assumes your restart file has velocities.

To run, just give it the paths to your prmtop, restart data and desired output data as below.

main('test_data/parm.prmtop','test_data/','test_data/outs/')

Functions🔗

convert_geometries(prmpath: str, rstpath: str, outpath: str) 🔗

This function creates a directory structure and output to convert a set of geometries from Amber prmtop/rst format to FMS90 Geometry.dat format.

Parameters:

Name Type Description Default
prmpath str

the path to the prmtop file

required
rstpath str

the path to the restart file

required
outpath str

the path to the output files

required

Returns:

Type Description

All converted geometries in the path outpath/filename.

Source code in scripts/convert_to_geomdat.py
def convert_geometries(
    prmpath: str,
    rstpath: str,
    outpath: str,
):

    """This function creates a directory structure and output to convert a set of geometries from Amber prmtop/rst format to FMS90 Geometry.dat format.

    Params:
      prmpath: the path to the prmtop file
      rstpath: the path to the restart file
      outpath: the path to the output files

    Returns:
      All converted geometries in the path outpath/filename.
    """

    findme = "*/*.rst*"

    Prmpath = Path(prmpath)
    Rstpath = Path(rstpath)
    Outpath = Path(outpath)

    # Get initial information from parmtop
    names, masses = read_prmtop(Prmpath)

    # Find all restart files in test data
    filenames = list(Rstpath.glob(findme))

    # Create new dats
    for filename in filenames:
        coords, vel, numats = read_rst(filename)
        output_geom(names, masses, coords, vel, numats, Outpath, filename)

output_geom(namelist: <built-in function array>, atommass: <built-in function array>, atomcoord: <built-in function array>, atomvel: <built-in function array>, atomnum: int, outpath: str, filename: str) 🔗

This function takes output from read_rst and read_prmtop as input and creates a Geometry.dat file, with appropriate numerical conversions to Bohr coordinates and a.u. time, to be used for AIMS simulations.

Parameters:

Name Type Description Default
namelist <built-in function array>

the vector of atom names as strings

required
atommass <built-in function array>

the vector of atomic masses as floats

required
atomcoord <built-in function array>

the 2d N by 3 array of atomic xyz coordinates

required
atomvel <built-in function array>

the 2d N by 3 array of atomic velocities

required
atomnum int

the total number of atoms in the system

required
outpath str

the directory to put all output files/directories

required
filename str

the name of the specific restart file to be converted

required

Returns:

Type Description

A Geometry.dat file in the path outpath/filename-Geometry.dat

Source code in scripts/convert_to_geomdat.py
def output_geom(
    namelist: np.array,
    atommass: np.array,
    atomcoord: np.array,
    atomvel: np.array,
    atomnum: int,
    outpath: str,
    filename: str,
):

    """This function takes output from read_rst and read_prmtop as input and creates a Geometry.dat file, with appropriate numerical conversions to Bohr coordinates and a.u. time, to be used for AIMS simulations.

    Params:
      namelist: the vector of atom names as strings
      atommass: the vector of atomic masses as floats
      atomcoord: the 2d N by 3 array of atomic xyz coordinates
      atomvel: the 2d N by 3 array of atomic velocities
      atomnum: the total number of atoms in the system
      outpath: the directory to put all output files/directories
      filename: the name of the specific restart file to be converted

    Returns:
      A Geometry.dat file in the path outpath/filename-Geometry.dat

    """

    # Initialize new output directory if it doesn't exist

    if not outpath.exists():
        Path.mkdir(outpath)

    # Create directory structure from original files

    newdir = filename.parent
    finalpath = outpath.joinpath(newdir.relative_to(newdir.parent))
    Path.mkdir(finalpath, parents=True, exist_ok=True)

    with open(finalpath.joinpath("Geometry.dat"), "w+") as f:
        f.write("UNITS=BOHR\n")
        f.write(str(atomnum) + "\n")
        for j in range(0, atomnum):
            if namelist[j] == "Cl-":
                f.write(
                    "{0:<10s}".format(namelist[j][0:2])
                    + "".join(
                        "{0:>22.16f}".format(x / 0.52917724924) for x in atomcoord[j]
                    )
                    + "\n"
                )
            else:
                f.write(
                    "{0:<10s}".format(namelist[j][0])
                    + "".join(
                        "{0:>22.16f}".format(x / 0.52917724924) for x in atomcoord[j]
                    )
                    + "\n"
                )
        f.write("# momenta\n")
        for j in range(0, atomnum):
            for x in atomvel[j]:
                p = x * atommass[j] * 9.3499611711905144e-04 * 1.8228884855409500e03
                f.write(
                    "{0:>22.16f}".format(p)
                )  # time unit in rst7: 1/20.455 ps   length unit in rst7: angstrom    mass unit in amber: atomic mass unit
            f.write("\n")

read_prmtop(prmpath: str) 🔗

This function parses a prmtop to extract atomic types/names and masses.

Parameters:

Name Type Description Default
prmpath str

the path to the prmtop file

required

Returns:

Type Description
namelist

the atom names as a vector of strings. atommass: atomic masses as a vector of floats.

Source code in scripts/convert_to_geomdat.py
def read_prmtop(prmpath: str):

    """This function parses a prmtop to extract atomic types/names and masses.

    Params:
      prmpath: the path to the prmtop file

    Returns:
       namelist: the atom names as a vector of strings.
       atommass: atomic masses as a vector of floats.

    """

    namelist = []
    atommass = []
    with open(prmpath) as f:
        for line in f:
            if "ATOM_NAME" in line:
                line = next(f)
                for line in f:
                    if "CHARGE" in line:
                        break
                    else:
                        line = line.strip("\n")
                        for c in range(0, len(line), 4):
                            namelist = namelist + [line[c : c + 3]]
            if "MASS" in line:
                line = next(f)
                for line in f:
                    if "ATOM_TYPE_INDEX" in line:
                        break
                    else:
                        a = line.split()
                        for x in a:
                            atommass.append(float(x))
    return namelist, atommass

read_rst(rstpath: str) 🔗

This function parses a restart file to extract coordinates and velocities.

Parameters:

Name Type Description Default
rstpath str

the path to the restart file(s) to be processed. Can handle instances in subdirectories.

required

Returns:

Type Description
atomcoord

the xyz coordinates as an N by 3 array. atomvel: the velocities as an N by 3 array.

Source code in scripts/convert_to_geomdat.py
def read_rst(rstpath: str):

    """This function parses a restart file to extract coordinates and velocities.

    Params:
      rstpath: the path to the restart file(s) to be processed. Can handle instances in subdirectories.

    Returns:
      atomcoord: the xyz coordinates as an N by 3 array.
      atomvel: the velocities as an N by 3 array.

    """

    atomcoord = np.array([])
    atomvel = np.array([])
    atomnum = 0
    with rstpath.open() as f:
        f.readline()
        line = f.readline()
        atomnum = int(line.split()[0])
        atomcoord = np.zeros((atomnum, 3))
        atomvel = np.zeros((atomnum, 3))
        for j in range(0, int((atomnum + 1) / 2)):
            line = f.readline()
            a = line.split()
            atomcoord[j * 2] = np.array([float(a[0]), float(a[1]), float(a[2])])
            if len(a) == 6:
                atomcoord[j * 2 + 1] = np.array([float(a[3]), float(a[4]), float(a[5])])
        for j in range(0, int((atomnum + 1) / 2)):
            line = f.readline()
            a = line.split()
            atomvel[j * 2] = np.array([float(a[0]), float(a[1]), float(a[2])])
            if len(a) == 6:
                atomvel[j * 2 + 1] = np.array([float(a[3]), float(a[4]), float(a[5])])
    return atomcoord, atomvel, atomnum