Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add load_final_frame to analysis module to facilitate MC-MD workflows #97

Open
rsdefever opened this issue May 5, 2021 · 1 comment

Comments

@rsdefever
Copy link
Collaborator

rsdefever commented May 5, 2021

Is your feature request related to a problem? Please describe.
Right now the user would need to write a custom function to read the final frame of a Cassandra trajectory.

Describe the solution you'd like
An API that looks like:

import mosdef_cassandra as mc
...
mc.load_final_frame(system, trajectory_filename)
mc.load_frame(system, trajectory_filename, frame_number)  # alternative

** Open questions **
What data structure should this return? mbuild.Compound? parmed.Structure? gmso.Topology?
The Compound is easiest to implement, but if we are going to require the system object for this function call, then perhaps we can return the full system/topology?

Describe alternatives you've considered
Another option is to create a to_mdtraj function -- but this requires adding mdtraj as a dependency.

@rsdefever
Copy link
Collaborator Author

Example code:

def load_final_xyz_frame(fname):
    """Return the final frame of a Cassandra .xyz file as an mbuild.Compound

    Assumes there is a .H file with the same name. E.g., if the .xyz file
    is 'equil.out.xyz', there should also be an 'equil.out.H' containing
    the box information.

    Arguments
    ---------
    fname : str
        path to the xyz file

	Returns
	-------
    mbuild.Compound
        the mbuild.Compound for the final frame
    """
    if not isinstance(fname, str):
        raise TypeError("'fname' must be a string")
    if fname[-4:] == ".xyz":
        fname = fname[:-4]

    data = []
    with open(fname + ".xyz") as f:
        for line in f:
            data.append(line.strip().split())

    for iline, line in enumerate(data):
        if len(line) > 0:
            if line[0] == "MC_STEP:":
                natom_line = iline - 1

    final_frame = data[natom_line + 2 :]
    natoms = int(data[natom_line][0])
    with open(fname + "-final.xyz", "w") as f:
        f.write(f"{natoms}\nAtoms\n")
        for coord in final_frame:
            f.write(
                "{}\t{}\t{}\t{}\n".format(
                    coord[0], coord[1], coord[2], coord[3],
                )
            )
    data = []
    with open(fname + ".H") as f:
        for line in f:
            data.append(line.strip().split())

    nspecies = int(data[-1][0])
    box_matrix = np.asarray(
        data[-(nspecies + 5) : -(nspecies + 2)], dtype=np.float32
    )
    assert box_matrix.shape == (3, 3)
    if np.count_nonzero(box_matrix - np.diag(np.diagonal(box_matrix))) > 0:
        raise ValueError("Only orthogonal boxes are currently supported")

    # If all is well load in the final frame
    frame = mbuild.load(fname + "-final.xyz")
    # mbuild.Compounds use nanometers!
    frame.periodicity = np.diagonal(box_matrix / 10.0)

    return frame

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant