About this package

A python package containing an MDAnalysis class to calculate charge, dipolar and quadrupolar density profiles from simulation trajectories. Also includes scripts (in the bin folder) to use this class from the command line and a class for calculating per molecule electrostatic properties.

If you use this package, please link to this page, cite our paper (Chapman and Bresme 2022, unpublished) and cite the papers of the MDAnalaysis authors

Requirements and Dependencies

Python Packages:

• setuptools
• MDAnalysis(https://www.mdanalysis.org/)
• numpy and matplotlib
• (optional) pip

How to install

These instructions should work on Linux and MacOS. You may need to do additional steps to add the scripts to your path on Windows.

To install, download the repo and change directory to it. Next either run:

• python setup.py install or
• pip install ./

The latter might be preferred as it is easier to uninstall (with pip uninstall multipole) and will try to automatically download dependancies, if available on pip. When using the first option you may want to use an enviroment manager, such as the one in conda.

Installing with previous instructions will add the script calc_multipoles_dens to your python installation's bin folder (e.g. for anaconda, ~/anaconda3/bin/). This is typically in your path, so you can use calc_multipoles_dens directly from the command line.

Running

The command line utility calc_multipole_dens has only one mandatory option, a topology file. This however will only give you the charge densities for the frame represented by that topology with the default settings. It is very important to change the settings for your needs.

The help message for calc_multipole_dens is as follows:

usage: calc_multipoles_dens [-h] [-f Trajectory] [-o Output] [-v] [-m]
                            [--water_model {tip4p05,spce,tip3p,opc}]
                            [-b BEGIN] [-e END] [-s STEP] [-M M_NAME]
                            [-H H_NAME [H_NAME ...]] [-d] [-w BIN_WIDTH] [-c]
                            [-u] [--types_or_names {None,type,name}]
                            topfile

Calculate the bins of the charge, dipole and quadrapole.Currently only
implemented for water.This can be done for just a single frame, or an entire
trajectory

positional arguments:
  topfile               The topology file name, .data file, for

optional arguments:
  -h, --help            show this help message and exit
  -f Trajectory, --trajfile Trajectory
  -o Output, --outfile Output
                        The output file name
  -v, --verbose         Increase verbosity
  -m, --inmem           Load the enture Trajectory into memory at once. May be
                        faster once loaded, but requires a lot of memory.
  --water_model {tip4p05,spce,tip3p,opc}
                        calculate the mutipole moments for the specified water
                        model. Applies some pre-processing such as assigning
                        charges if not already.TODO: Don't assign charges
                        Future versions will make None default, so non water
                        systems can be used.
  -b BEGIN, --begin BEGIN
                        The first frame to use
  -e END, --end END     The last frame to use. -1 means the last frame
  -s STEP, --step STEP  The step size to use
  -M M_NAME, --M_name M_NAME
                        the atom name or type to use for the centre of the
                        molecule for binning and multipole moment calculation.
                        Should use O for tip3p/spc styles and dummy atom for
                        tip4p styles, or Oxygen if calculating dummy atom on
                        the fly.
  -H H_NAME [H_NAME ...], --H_name H_NAME [H_NAME ...]
  -d, --calculate_dummy
                        Calculate the position of the dummy atom on the fly,
                        from the positions of the oxygen and hydrogen atoms.
  -w BIN_WIDTH, --bin_width BIN_WIDTH
                        The target width of the bins to use, in Angstrom.
                        Overridden if number of bins is set instead
  -c, --coord_centre    if used, binned coordinate is centre of the bin, else
                        it is the left edge
  -u, --unwrap          Unwrap the coordinates
  --types_or_names {None,type,name}
                        Specify whether hydrogen/centre choice is in reference
                        to the atom type or name. By default tries to work it
                        out from the topology file.

The option -f is where you specify the trajectory file you wish to use. It is important that you set the water model to match that of the system you are using, with the --water_model option. Currently only the TIP4P/2005, SPC/E, TIP3P and OPC models work.

You need to set -M and -H to the name or type of the negatively charged atom (oxygen or the dummy atom) and the hydrogen atoms, respectively. Whether name or type is used depends on the --types_or_names option.

Because this package is based on MD analysis, it should accept a wide range of trajectory and topology file formats, however, not all have been tested to work with this script. The topology and trajcetory formats MD analysis can accept are listed here.

Understanding the results

Running the script will output a file (that can be read with np.load_txt, for example) containing the charge densites, the multipole densities, the molecule density and the orientations of the molecules.

To calculate the the electostatic field and potentials, and dipolar and quadrupolar contributions, you need to integrate the charge/multipolar densities. The package thermotar has some tools for doing this.

Limitations

• Currently these scripts work with only a few water models.
• ~~If your water model uses dummy atoms, but the trajectories do not include the positions of the dummy atoms (such as the case with the output from LAMMPS), you must add these back in. This script does not yet do this whilst reading the trajectory, however you can generate a trajectory that re-adds the atoms using our dummify package. ~~
• The location of the origin for calculating the dipole and quadrupole of each molecule must be on the negatively charged atom (real or dummy). The choice in orign has no effect on the dipole moments, but can affect the quadrupole contribution. For water the location of the negative charge is a good choice.
• Electrostatic fields/potentials must be calculated from the average density profiles in post processing, no support for calculating for each frame yet.
• No sub-averaging for estimating errors in one trajectory is performed. Use independant repeats and average and calculate the standard error of the potential profiles from these repeats.
• Binning is only avalible in 1D currently.
• Charges that change in time are unsupported.

Changes From 0.0.1 to 0.1

• Can calculate the position of the dummy atoms on the fly using the -d option.
• Performs unwrapping of water molecules on the fly (needed for multipole and dummy calculation). Specify with -u.
• Supports the presence of charges not due to water molecules. Binned in the total charge and independently.

• Now requires numba for fast calculations of the quadrupole, dipole moments, box unwrapping and dummy atom position.
  • Due to JIT compilation, This may be slower for short trajectories/single frames. In the future this may be optional.
• Now requires fast_histogram for fast calculation of the charge binning. This was chosen because the numba implementation of np.histogram does not support weights, which are essential for binning.

Warning about DCD Trajectories

If using MDAnalysis version 2.4.X for X < 3, there is a bug when reading DCD large trajectories that causes crashing. This is fixed in newer versions, so make sure to use the latest version or use 2.3 or earlier.

TODO

• Unwrap trajectories using oxygen hydrogen distance, or MD ANAlysis residues.
• Add the option to only calculate the charge profiles.
• Add calculations of the potential and field at the end of the simulation.
• Update the README.md for the command line arguments.
• Calculate charge densities in e/AA^3. Add option for old behaviour.