hougroup / mdgo

A codebase for classical molecular dynamics (MD) simulation setup and results analysis.

An all-in-one code base for the classical molecualr dynamics (MD) simulation setup and results analysis.

1. Installation

1.1 Installing from PyPI

To install the latest release version of mdgo:

pip install mdgo

1.2 Installing from source code

Mdgo requires numpy, pandas, matplotlib, scipy, tqdm, statsmodels, pymatgen>=2022.0.9, pubchempy, selenium, MDAnalysis (>=2.0.0) and their dependencies.

Getting Source Code

If not available already, use the following steps.

1. Install git, if not already packaged with your system.

2. Download the mdgo source code using the command:

   https://github.com/HouGroup/mdgo.git

Installation

1. Navigate to mdgo root directory:

   cd mdgo

2. Install the code, using the command:

   pip install .

3. The latest version MDAnalysis (>=2.0.0) is recommended, you may download the source code of the latest MDAnalysis from github and install using pip to replace an existing version.

Installation in development mode

1. Navigate to mdgo root directory:

   cd mdgo

2. Install the code in "editable" mode, using the command::

   pip install -e .

3. The latest version MDAnalysis (>=2.0.0) is recommended, you may download the source code of the latest MDAnalysis from github and install using pip to replace an existing version.

2. Features

1. Retrieving compound structure and information from PubChem
   • Supported searching text:
     • cid, name, smiles, inchi, inchikey or formula
   • Supported output format:
     • smiles code, PDB, XML, ASNT/B, JSON, SDF, CSV, PNG, TXT
2. Retrieving water and ion models
   • Supported water models:
     • SCP, SPC/E, TIP3P_EW, TIP4P_EW, TIP4P_2005
   • Supported ion models:
     • alkali, ammonium, and halide monovalent ions by Jensen and Jorgensen
     • alkali and halide monovalent ions by Joung and Cheatham
     • alkali and alkaline-earth metal cations by Åqvist
3. Write OPLS-AA forcefield file from LigParGen
   • Supported input format:
     • mol/pdb
     • SMILES code
   • Supported output format:
     • LAMMPS(.lmp)
     • GROMACS(.gro, .itp)
4. Write OPLS-AA forcefield file from Maestro
   • Supported input format:
     • Any format that Maestro support
   • Supported output format:
     • LAMMPS(.lmp)
     • Others pending...
5. Packmol wrapper
   • Supported input format:
     • xyz
     • Others pending...
6. Basic simulation properties
   • Initial box dimension
   • Equilibrium box dimension
   • Salt concentration
7. Conductivity analysis
   • Green--Kubo conductivity
   • Nernst--Einstein conductivity
8. Coordination analysis
   • The distribution of the coordination number of single species
   • The integral of radial distribution function (The average coordination numbers of multiple species)
   • Solvation structure write out
   • Population of solvent separated ion pairs (SSIP), contact ion pairs (CIP), and aggregates (AGG)
   • The trajectory (distance) of cation and coordinating species as a function of time
   • The hopping frequency of cation between binding sites
   • The distribution heat map of cation around binding sites
   • The averaged nearest neighbor distance of a species
9. Diffusion analysis
   • The mean square displacement of all species
   • The mean square displacement of coordinated species and uncoordinated species, separately
   • Self-diffusion coefficients
10. Residence time analysis
    • The residence time of all species

3. Citation

When using mdgo in published work, please cite the following paper:

• Hou, T.; Fong, K. D.; Wang, J.; Persson, K. A. The solvation structure, transport properties and reduction behavior of carbonate-based electrolytes of lithium-ion batteries. Chem. Sci. 2021, 12, 14740-14751. [doi]

mdgo's Issues

Conversion of LJ parameters for ions

Hi @htz1992213 , I'm trying to understand how to convert literature Lennard Jones parameters into LAMMPS format, and I'm concerned that the ones for ions stored in /data/ions may not have been converted correctly.

For example, consider Na+ from the Joung-Cheatham paper here

Here is what we have in mdgo for Na+ for the SPCE model:
https://github.com/htz1992213/mdgo/blob/c843a9643a67992cf5069e7be668b194e86ef889/mdgo/data/ion/joung_cheatham/spce/Na%2B.lmp#L16-L18

Here is the relevant excerpt from the paper:

The LAMMPS docs say that

Note that sigma is defined in the LJ formula as the zero-crossing distance for the potential, not as the energy minimum at 2^(1/6) sigma

For whatever reason, it's conventional to report sigma as r_min/2 in many papers I've seen. The potential should equal zero when sigma = r. So if I understand correctly, to convert rmin/2 to the sigma that LAMMPS wants, you should do

(r_min / 2 * 2) / 2^(1/6)

If I perform this conversion on the Na+ value from the paper, I get 2.159, vs. 2.238 in mdgo. Am I missing something? Also, would there be benefit to storing more decimal places for the epsilon parameter?

This is incredibly tricky and I'm really surprised that this subtle point about how LAMMPS defines LJ parameters is not documented more clearly (or that the research community doesn't just report sigma instead of r_min/2). Please let me know what you think, because as I add new parameters I want to make sure I'm converting them correctly.

CombinedData vs. LammpsData.from_file causes MDRun.init to fail

I am trying to instantiate an MDRun from a topology file md.data (attached) and a trajectory in dcd format. Per the MDRun docstring, the first argument is supposed to be a LammpsData object. So, I load the object from my file:

data = LammpsData.from_file('md.data')

and then pass it to MDRun

tma_cl = MdRun(data,
             traj,
             traj_unwrapped,
             0, 1, "TMA-Cl",
             )

where traj and traj_unwrapped are wrapped and unwrapped dcd files, respectively. The init fails because of this check for CombinedData

-----------------------------------------------------------------------
AssertionError                         Traceback (most recent call last)
<ipython-input-20-feee45d507c3> in <module>
      2 from pymatgen.io.lammps.data import CombinedData, LammpsData
      3 # load lammps inputs
----> 4 tma_cl = MdRun(data,
      5              traj, traj_unwrapped,
      6              0, 1, "TMA-Cl",

~/miniconda3/envs/md/code/mdgo/mdgo/core.py in __init__(self, lammps_data, wrapped_run, unwrapped_run, nvt_start, time_step, name, select_dict, res_dict, cation_name, anion_name, cation_charge, anion_charge, temperature, cond)
    105         self.res_dict = res_dict
    106         if self.select_dict is None and self.res_dict is None:
--> 107             self.res_dict = res_dict_from_lammpsdata(self.data)
    108         if self.res_dict is None:
    109             self.res_dict = res_dict_from_select_dict(self.wrapped_run, self.select_dict)

~/miniconda3/envs/md/code/mdgo/mdgo/util.py in res_dict_from_lammpsdata(lammps_data)
    439         dict: A dictionary of resnames.
    440     """
--> 441     assert isinstance(lammps_data, CombinedData)
    442     idx = 1
    443     res_dict = dict()

If I disable the assert statement, I instead get

------------------------------------------------------------------------
AttributeError                         Traceback (most recent call last)
<ipython-input-21-1e8e201524a1> in <module>
      2 from pymatgen.io.lammps.data import CombinedData, LammpsData
      3 # load lammps inputs
----> 4 tma_cl = MdRun(data,
      5              traj,
      6              traj_unwrapped,

~/miniconda3/envs/md/code/mdgo/mdgo/core.py in __init__(self, lammps_data, wrapped_run, unwrapped_run, nvt_start, time_step, name, select_dict, res_dict, cation_name, anion_name, cation_charge, anion_charge, temperature, cond)
    105         self.res_dict = res_dict
    106         if self.select_dict is None and self.res_dict is None:
--> 107             self.res_dict = res_dict_from_lammpsdata(self.data)
    108         if self.res_dict is None:
    109             self.res_dict = res_dict_from_select_dict(self.wrapped_run, self.select_dict)

~/miniconda3/envs/md/code/mdgo/mdgo/util.py in res_dict_from_lammpsdata(lammps_data)
    457                     res_dict[n] = "same mass as resid " + str(start + i)
    458     else:
--> 459         for name, num in zip(lammps_data.names, lammps_data.nums):
    460             start = idx
    461             idx += num

AttributeError: 'LammpsData' object has no attribute 'names'

If I load md.data as a CombinedData and not a LammpsData

data = LammpsData.from_file('md.data')

I get

------------------------------------------------------------------------
TypeError                              Traceback (most recent call last)
<ipython-input-22-f46b835926f9> in <module>
----> 1 data = CombinedData.from_file(pipeline_dir / '2021-05_prelim_lammps_TMA-Cl' / 'md.data')

~/miniconda3/envs/md/code/pymatgen/pymatgen/io/lammps/data.py in from_file(cls, filename, atom_style, sort_id)
    773         items["atom_style"] = atom_style
    774         items["box"] = box
--> 775         return cls(**items)
    776 
    777     @classmethod

TypeError: __init__() got an unexpected keyword argument 'masses'

I think this is because CombinedData inherits its from_file method from LammpsData, and LammpsData.from_file does not accept the masses kwarg. What's confusing to me though is that if I load md.data as a LammpsData object, data.masses is populated correctly.

What is the correct way to fix this? Is this a problem with CombinedData.from_file() or with MDRun.init(), or am I misunderstanding how to pass data into MDRun?

Thanks for any assistance!

md.data.txt

Feature request: `step` kwarg for residence time related functions

Most core methods in MDAnalysis accept the kwargs start, stop, and step to specify the starting and ending frames of the trajectory to analyze, and how many frames to skip (i.e., only analyze every nth frame where n=step).

In mdgo, the residence time functions like get_neighbor_corr accept start and stop but not step. For RTD in particular, one usually wants to start with the first frame and analyze a long time into the simulation, which can take a very long time if every frame is considered. It would be great to be able to pass step to this method so that a preliminary analysis can be done more quickly.

There may be technical reasons why skipping frames is not advised and/or makes the analysis invalid. But unless that's the case, I think this would be a nice feature to add. I'm happy to help but may not have time to work on it for a while, so I thought I'd open this issue as a placeholder.

python support < 3.8

Hey @rkingsbury
I found an import error for Final type when using python version < 3.8. I have removed that Final typing for a quick fix. Any idea of how that could be better addressed?
https://github.com/htz1992213/mdgo/blob/06b50b28d805d795fe945c48e7d6faa01c02506c/mdgo/util.py#L13
https://github.com/htz1992213/mdgo/blob/06b50b28d805d795fe945c48e7d6faa01c02506c/mdgo/forcefield.py#L47