Construction of Multi-component Systems

Building basic configurations of multicomponent and/or multilayered systems normally proceed by initially placing the primary solute molecules, e.g. a protein or a crystallographic mineral slab, into a pre-defined simulation box using functions such as place_atom(), center_atom(), translate_atom(), wrap_atom(), and so on.

Ions can be added into a simulation box either randomly in specified regions (but without atomic overlap of existing molecules), on specified planes (in 3D), or with exponentially decaying concentrations from a surface with the function ionize_atom(). Other similar functions to add new or copy existing solutes also exist.

In terms of solvation, the function solvate_atom() can solvate an entire simulation box or limited regions with generic three-, four-, or five-site water models, such as SPC/E or TIP3P/TIP4P/TIP5P, either as liquid water or hexagonal ice. Solvation with custom solvents using unwrapped solvent boxes is also supported, as is scaling the solvent density to optimize solvent molecule packing. Furthermore, one can solvate a solute by a certain solvation thickness (Ex. around a centered protein).

Visualization

For visualization one can make use of the function vmd() if the VMD software is installed separately (Humphrey et al. Reference Humphrey, Dalke and Schulten1996), and the PATH2VMD() function is properly set. Alternatively, the basic plot_atom() or the plot_density_atom() functions rapidly plot a full simulation box with thousands of atoms (see Fig. 2), with or without the density profiles along the x,y,z-dimensions. Analogously, and although it is significantly slower the show_atom() and the show_density_atom() functions can be used to render glossy atoms with cylindrical bonds or VdW spheres.

Fig. 2 Example image of a 136,000 particle system generated with the plot_density_atom() function, showing scaled density profiles of a multilayered montmorillonite particle (black profiles), Na⁺ counter-ions (blue profiles) hydrated by three pseudo-layers of water (O_water in red, H_water in gray profiles)

Bonding and Structural Analysis

Most functions, such as dist_matrix_atom() and bond_angle_atom(), take PBC into account, which allows for generation of molecular topologies with bonds, angles and basic dihedrals across the PBC. Many functions also support triclinic boxes using the ‘tilt vectors’ xy, xz, yz as defined in Gromacs and LAMMPS. Several of the functions use two types of neighbor cutoffs to find interacting atoms, one shorter cutoff for bonded H’s (default 1.25 Å) and one larger cutoff for all other interactions (default 2.25 Å). Furthermore, a general neighbor-distance analysis can be performed: (1) by comparison with the Revised Shannon radii (Shannon Reference Shannon1976; van Horn Reference van Horn2001); (2) by calculating the bond-valence values using the semi-empirical Bond Valence Method (Brown Reference Brown2009, Reference Brown and Mingos2013, Reference Brown2016) or; (3) by calculating basic XRD profiles. Here, analysis 2 is performed most easily by calling the analyze_atom() function, which first finds all atom bonds and angles to the nearest neighbors. Secondly, it makes educated guesses based on the chemical name and apparent coordination number, the oxidation state, electron configuration, and the ideal crystalline, ionic, and vdW-radii of the atoms. This function is, hence, useful for sanity-checks of either the initial input structures or to scrutinize one simulation result. As the Bond Valence Method also calculates the Global Instability Index (G _II), it can further be used in structural data-mining studies to check for structural stress in molecules. As for analysis 3 and XRD analysis, a basic XRD function called XRD_atom() for single-crystal UCs is provided, using 11-coefficient atomic scattering factors (Waasmaier & Kirfel Reference Waasmaier and Kirfel1995).

Trajectory Analysis

Although not the primary purpose of the atom library, support for trajectory reading (and writing) of binary trajectories utilizing the third-party MatDCD or the mxdrfile packages for the .dcd and .xtc|.trr trajectory formats, respectively, is provided (Gullingsrud Reference Gullingsrud2000; Kapla & Lindén Reference Kapla and Lindén2018). For reading .xtc trajectories, one could also use the Gro2mat package (Dien et al. Reference Dien, Deane and Knapp2014). In addition to these third-party binary trajectory packages, text-format trajectories in the .pdb|.xyz|.gro formats are also natively supported by the atom library, where one special version written specifically for the MC code RASPA2 named import_mc_pdb_traj() can handle trajectories with non-constant numbers of atoms.

In order to import a trajectory into MATLAB as shown below, both a coordinate file and a trajectory file is needed.

$[atom,traj]=import\_traj\text{(}...coordfile,trajfile)$

Apart from the atom struct taken from the coordinate file, this command results in a data matrix called traj, of size 3N×frames, where N represents the number of atoms in the trajectory.

Case Study: A Hydrated Na⁺-montmorillonite Nanopore System with CLAYFF atomtypes

This example demonstrates the main steps needed to setup a simulation box for MC or MD simulations of a system representing the water–mineral interface to a defect mineral particle called MMT (as in Montmorillonite). MMT will describe a negatively charged smectite layer, with an approximate thickness of 1 nm and being infinite in the lateral plane, due to the fact that PBC is used in MC or MD simulations and because the layer in this example will have no edges in order to avoid complexity.

Preparation of the Mineral Lattice

The negative charge of MMT mainly originates from isomorphic substitution of octahedrally coordinated Al(III) with Mg(II), which creates a charge defect which is smeared over the nearest coordinating O(-II) in the lattice. In most cases when performing molecular modeling of inorganic or mineral particles, initial UC data are available from either X-ray diffraction data or quantum mechanical calculations. For MMT, due to the isomorphic substitution, no precise UC data exist. Instead, one can use isostructural and defect-free UC data from pyrophyllite (Lee & Guggenheim Reference Lee and Guggenheim1981), and then replicate it to a mineral layer and subject it to isomorphic substitution. Hence, initially, a pyrophyllite UC structure must be imported:

$Pyrox1=import\_atom(\text{'}Pyrophyllite.pdb\text{'})$

In many cases when UC information is taken from XRD data, the number and positions of H atoms are somewhat unreliable. Hence, in this example both initial removal and addition of H atoms to the clay lattice are demonstrated. First, all atoms with names starting with ‘H’ are removed.

$Pyrox1\left(strncmpi([Pyrox1.type],\text{'}H\text{'},1)\right)=\left[\right]$

Then the remaining mineral lattice is analyzed for the resulting structure using the function analyze_atom().

$properties=analyze\_atom(Pyrox1,Box\_\dim )$

This function will output a separate properties struct with bonding information, as well as reveal which atoms need healing by outputting a heal_ind variable. To protonate those individual atoms, one would issue:

$H\_atom=protonate\_atom\text{(}...Pyrox1,Box\_\dim ,heal\_ind\text{)}$

To simplify this example one could convert the triclinic UC to an orthogonal UC:

$Pyrox1=orto\_atom\text{(}...Pyrox1H\_atom,Box\_\dim \text{)}$

Next, the UC must be replicated into a layer, for instance by 6×4 times in the X and Y-directions, respectively:

$Pyrox6x4=replicate\_atom\text{(}...Pyrox1,orto\_Box\_\dim ,\left[641\right]\text{)}$

From this pyrophyllite layer, a MMT layer can be generated by performing isomorphic substitution by replacing Al lattice atoms with Mg in two-thirds of all the UCs (which is typical for natural MMT clay minerals). The last numerical argument (in Ångström units) sets the nearest Mg–Mg neighbor distance, as charge defects are unlikely to be close to each other.

$MMT=substitute\_atom\text{(}Pyrox6x4,...Box\_\dim ,6\ast 4\ast 2/3,\text{'}Al\text{'},\text{'}Mg\text{'},5.5\text{)}$

For consistency, one can at this point change the attribute MMT.resname from ‘PYR’ to ‘MMT’:

$\left[MMT.resname\right]=deal\left(\left(\text{'}MMT\text{'}\}\right)\right)$

Here, the deal() is a native MATLAB function. Next, the names of all atoms (i.e. the atomtypes) in MMT can be set to a specific forcefield, such as the original CLAYFF (Cygan et al. Reference Cygan, Liang and Kalinichev2004).

$MMT=clayff\_2004\_atom(MMT,Box\_\dim ,\text{'}clayff\text{'})$

This function will also output the layer net charge (–16 from the isomorphic substitution), which must be counterbalanced by counter-cations such as Na⁺. In order to accommodate these ions and later water molecules, the dimensions of the definitive simulation box should be defined, from the lateral dimensions of the MMT layer along X/Y (already set in the existing Box_dim variable) and an arbitrary height in Z, here set to 40 Å:

$Box\_\dim \left(3\right)=40;\%inAngstrom$

Adding Counterions and Solvating with Water

Ions can be added to an existing simulation box in several ways using the ionize_atom() function. In this example, an Ions struct is created by adding ions close to the mineral surface (but >3 Å) of the wrapped MMT:

$Ions=ionize\_atom(...\text{'}Na\text{'},\text{'}Na\text{'},Box\_\dim ,16,3,...wrap\_atom(MMT,Box\_\dim ),\text{'}surface\text{'})$

Note that wrapping the MMT layer into the box before adding the ions is important in order to avoid atomic overlap. Next, a Water struct is created in order to solvate the simulation box. This is accomplished by placing 1000 SPC/E water molecules with a relative density of 1.05 no nearer than 2 Å to any solute:

$Water=solvate\_atom(Box\_\dim ,1.05,...2,1000,update\_atom\left(\left(MMT,Ions\right)\right))$

Finally, all components are merged together using update_atom() into single System struct, which updates all atom- and molecular-ID indexes, and an ouput .pdb file is written. Note that because the different System struct components are stored in separate struct variables, the order of the components can be changed arbitrarily.

$System=update\_atom\left(\left\{MMTWaterIons\right\}\right)write\_atom\_pdb(System,...Box\_\dim ,\text{'}System.pdb\text{'})$

The following lines show example commands used to plot the final System struct (see Fig. 3), where, for instance, the optional last numerical in the plot_atom() function denotes the VdW radii scaling factor, whereas the 'vdw' in the show_atom() function sets the display style:

Fig. 3 The final hydrated and unwrapped Na⁺-MMT System generated by the plot_atom(). Note the legend which indicates the CLAYFF assigned atom types for the Water, Ions, and MMT struct variables. Alternatively, one could plot the System struct or the printed .pdb file with VMD instead.

$plot\_atom(System,Box\_\dim ,0.2)show\_atom(System,Box\_\dim ,\text{'}vdw\text{'})vmd(System,Box\_\dim )vmd(\text{'}System.pdb\text{'})$

Writing Molecular Topology Files

The atom library contains different flavors of both the CLAYFF and INTERFACE forcefields. The following examples show how to write basic .psf, .lmp, and .itp topology files according to the original CLAYFF to be used in the MD software like NAMD2, LAMMPS, and Gromacs.

$write\_atom\_psf(System,Box\_\dim ,\text{'}System\text{'})write\_atom\_lmp(System,Box\_\dim ,\text{'}System\text{'})write\_atom\_itp(MMT,Box\_\dim ,\text{'}MMT\text{'})$

Note that when generating molecular topology files, one must, in general, decide which components to include (i.e. the MMT, Ions, and Water) in the molecular topology, as different MD and MC software may use different approaches on how to include molecular topologies in the total system topology. For instance, in Gromacs, separate molecular topology files are used commonly (.itp files) for different components; hence, the Ions and the Water in the current example would preferably not be included in the same molecular topology as MMT.

On a separate note, because the original CLAYFF (nor INTERFACE) does not contain all different oxygen atom types needed to model all minerals or, for instance, custom clay edge models, the atom library can also be used to derive partial charges and molecular topologies for new oxygen atom types not part of the original CLAYFF forcefield. This is accomplished by equation 1, which in fact can be used to calculate the partial charge of any new or original oxygen atom type in CLAYFF (Tournassat et al. Reference Tournassat, Bourg, Holmboe, Sposito and Steefel2016; Lammers et al. Reference Lammers, Bourg, Okumura, Kolluri, Sposito and Machida2017).

(1) $Z_O=-2-(-\sum _n\frac{Z_n^f-Z_n^p}{C{N}_n})$

In this equation, the partial charge of a CLAYFF oxygen type (Z _O) is obtained from its formal charge (–2) minus the total charge needed to balance (hence the extra minus sign) the charge distributed over all the n nearest coordinating cationic centers (including H). This distributed charge represented by the second term on the right hand side is, in turn, calculated from the difference between the formal and partial charge $Z_n^f$ and $Z_n^p$ of each cation (hence the summation), divided by their respective number of coordinating oxygens, CN _n . For a more extensive explanation, see Lammers et al. (Reference Lammers, Bourg, Okumura, Kolluri, Sposito and Machida2017).

Analyzing the Lattice Structure

In order to study the structural integrity of molecules subjected to molecular simulations, computing the root mean square deviation of atomic positions (RMSD) from the positions in the input structure is an often-used method. This method is frequently used for large biomolecules; however, for crystalline inorganic lattice structures the semi-empirical Bond Valence Method can also be used. This method finds strained atoms in a lattice by relating the so-called bond valence values to the ideal semi-empirical bond distances (Brown Reference Brown2016), as well as by comparing the total bond valence sum to the ideal atomic valence. In this example, this method was applied by invoking the analyze_atom() function on the average MMT lattice structure (Fig. 3) obtained from equilibrated simulations performed in Gromacs in the NVT ensemble over 1 ns. Six different forcefield implementations were tested in total, namely: (1–2) two implementations of the INTERFACE forcefield (Heinz et al. Reference Heinz, Koerner, Anderson, Vaia and Farmer2005, Reference Heinz, Lin, Kishore Mishra and Emami2013), using the experimental or pre-defined bond distance (with the 5% scaling factor) and angle values, respectively; (3) the original unconstrained CLAYFF (Cygan et al. Reference Cygan, Liang and Kalinichev2004), having no defined metal–oxygen bonds or angles; (4) CLAYFF with no bonds but all M–O angle terms defined by the experimental values (as in 1); (5) a modified version of CLAYFF optimized for 1D-XRD modeling, using a larger basal O radius but smaller Si radius (Szczerba & Kalinichev Reference Szczerba and Kalinichev2016); and lastly (6) the same modified version of CLAYFF, but with all angle terms defined as in (4).

The resulting average bond distances (Fig. 4A) from the simulations generally agreed well among the various INTERFACE forcefields, apart from the Mg bond distance for CLAYFF, which was found to relax considerably compared to Al and the INTERFACE forcefields.

Fig. 4 Results from the structural analysis comparing six different implementations of the INTERFACE and CLAYFF forcefields (as explained in the text). (a) Average bond distances for the generic clay lattice atomtypes. (b) The corresponding bond valence sum, compared to the ideal atomic valence, shown as gray lines

Average bond distances reveal little, however, about the strain in lattices. The extended Mg–O bond distances found with the unconstrained CLAYFF but not with the constrained INTERFACE forcefields is reasonable considering the bond valence sum values. This is because a relaxed coordinating octahedron of the neighboring O atoms is needed in order to maintain a physically reasonable atomic valence (+2 for Mg, see Fig. 4B). Nevertheless, the Global instability index, G _II, which is a measure of the structural strain for the six different forcefield implementations (1–6) were 0.257, 0.230, 0.366, 0.233, 0.294, and 0.237, respectively. Overall this indicates that the implementations of the forcefields used are not perfect, as G _II values >0.2 from experimentally determined structures are found rarely (Brown Reference Brown2009), and hence indicate unstable lattices. Simulations with pyrophyllite (not shown) resulted in similar G _II values, which can be compared to the input pyrophyllite structure having a G _II value of 0.09. Although the INTERFACE 2013 implementation resulted in the lowest G _II value and, hence, the best representation of the internal atomic UC positions, the CLAYFF implementations with defined angles perform almost as well, and are better suited for 1D-XRD modeling due to the larger radius of the basal O atoms (Ferrage et al. Reference Ferrage, Sakharov, Michot, Delville, Bauer and Lanson2011; Szczerba & Kalinichev Reference Szczerba and Kalinichev2016). Still, the radii of the atom types in CLAYFF could likely be further optimized to better represent the bond distances within the actual UC.