Topology

The topology describes atomic and molecular properties as well as processes and reactions.

Global Properties

The following keywords control temperature, simulation box size etc., and must be placed outer-most in the input file.

temperature: 298.15  # system temperature (K)
geometry:
  type: cuboid       # Cuboidal simulation container
  length: [40,40,40] # cuboid dimensions (array or number)
mcloop:              # number of MC steps (macro × micro)
  macro: 5           # Number of outer MC steps
  micro: 100         # Number of inner MC steps; total = 5 × 100 = 500
random:              # seed for pseudo random number generator
  seed: fixed        # "fixed" (default) or "hardware" (non-deterministic)

Geometry

Below is a list of possible geometries, specified by type, for the simulation container, indicating if and in which directions periodic boundary conditions (PBC) are applied. Origin ($0,0,0$) is always placed in the geometric center of the simulation container. Particles are always kept inside the simulation container with an external potential that is zero if inside; infinity if outside.

geometry PBC Required keywords
type:
cuboid $x,y,z$ length (array or single float)
slit $x,y$ length (array or single float)
hexagonal $x,y$ radius (inscribed/inner), length (along z)
cylinder $z$ radius, length (along z)
sphere none radius

Simulation Steps

The variables macro and micro are positive integers and their product defines the total number simulations steps. In each step a random Monte Carlo move is drawn from a weighted distribution. For each macro step, all analysis methods are, if befitting, instructed to flush buffered data to disk and may also trigger terminal output. For this reason macro is typically set lower than micro.

Atom Properties

Atoms are the smallest possible particle entities with properties defined below.

atomlist Description
activity=0 Chemical activity for grand canonical MC [mol/l]
pactivity −log10 of chemical activity (will be converted to activity)
alphax=0 Excess polarizability (unit-less)
dp=0 Translational displacement parameter [Å]
dprot=0 Rotational displacement parameter [radians]
eps=0 Lennard-Jones/WCA energy parameter [kJ/mol]
mu=[0,0,0] Dipole moment vector [eÅ]
mulen=|mu| Dipole moment scalar [eÅ]
mw=1 Molecular weight [g/mol]
q=0 Valency or partial charge number [$e$]
r=0 Radius = sigma/2 [Å]
sigma=0 2r [Å] (overrides radius)
tension=0 Surface tension [kJ/mol/Å$^2$]
tfe=0 Transfer free energy [kJ/mol/Å$^2$/M]
psc Patchy sphero-cylinders properties (object)

A filename (.json) may be given instead of an atom definition to load from an external atom list. Atoms are loaded in the given order, and if it occurs more than once, the latest entry is used.

Example:

atomlist:
  - Na: {q:  1.0, sigma: 4, eps: 0.05, dp: 0.4}
  - Ow: {q: -0.8476, eps: 0.65, sigma: 3.165, mw: 16}
  - my-external-atomlist.json
  - ...

Anisotropic particles (experimental)

Faunus supports anisotropic particles such as dipoles, spherocylinders, and quadrupoles. While the implementation is pretty stable, it is not well optimized and should be considered experimental.

Patchy Sphero-cylinders

Sphero-cylinders are activated using the psc keyword. Currently only a single pair-potential is readily available, namely a cosine attraction for patch-patch interactions, mixed with Weeks-Chandler-Andersen for excluded volume.

atomlist:
  - mycigar:
      psc: {type: capped, length: 40, patch_angle: 80}
      eps: 2.4  # used for both wca and the cosine attraction
      sigma: 10 # cylinder and cap radius
      rc: 11.2  # cosine attraction switch distance
      wc: 6.0   # cosine attraction switch width

energy:
  - nonbonded:
      default:
        - coswca-psc:          # cosine attraction + WCA pair potential
            cos2: {mixing: LB} # used for interactions with patch
            wca: {mixing: LB}  # used for interactions with cylinder

Possible patch types are none, full and capped signifying no patch (○◻○); a patch that runs the full length including the end caps (●◼●); or a capped variant where the patch stops before the caps (○◼○). The coswca-psc potential can handle mixtures of PSCs and spherical particles. When a spherical particle interacts with the PSC, the closest distances to the patch and excluded volume are used to evaluate the attractive and repulsive energy.

psc Description
type=none none (○◻○), full (●◼●), capped (○◼○)
length=0 Length of cylinder (Å)
chiral_angle=0 Patch angle with respect to length axis (degrees)
patch_angle=0 Opening angle of patch (degrees), ◔ (degrees)
patch_angle_switch=0 (degrees)

Molecule Properties

A molecule is a collection of atoms, but need not be associated as real molecules. Two particular modes can be specified:

  1. If atomic=true the atoms in the molecule are unassociated and is typically used to define salt particles or other non-aggregated species. No structure is required, and the molecular center of mass (COM) is unspecified.
  2. If atomic=false the molecule resembles a real molecule and a structure or trajectory is required.

Properties of molecules and their default values:

moleculelist Description
activity=0 Chemical activity for grand canonical MC [mol/l]
atomic=false True if collection of atomic species, salt etc.
atoms=[] Array of atom names; required if atomic=true
bondlist List of internal bonds (harmonic, dihedrals etc.)
compressible=false If true, molecular internal coordinates are scaled upon volume moves
ensphere=false Radial rescale of positions to sphere w. radius of average radial distance from COM (stored in 1st atom which is a dummy)
excluded_neighbours=0 Generate an exclusionlist from the bonded interaction: Add all atom pairs which are excluded_neighbours or less bonds apart
exclusionlist List of internal atom pairs which nonbonded interactions are excluded
implicit=false Mark as implicit for reactive Monte Carlo schemes
insdir=[1,1,1] Insert directions are scaled by this
insoffset=[0,0,0] Shifts mass center after insertion
keeppos=false Keep original positions of structure
keepcharges=true Keep charges of structure (aam/pqr files) and traj even if mismatch with those in atomlist
rigid=false Set to true for rigid molecules. Affects energy evaluation.
rotate=true If false, the original structure will not be rotated upon insertion
structure Structure file or direct information; required if atomic=false
to_disk=false Save initial structure to {name}-initial.pqr; for molecular groups only
traj Read conformations from PQR trajectory. Cannot be used w. structure; see also keepcharges
trajweight One-column file with relative weights for each conformation. Must match frames in traj file.
trajcenter=false Move CM of conformations to the origin assuming whole molecules

Example:

moleculelist:
  - salt: {atoms: [Na,Cl], atomic: true}
  - water:
      structure: water.xyz
      bondlist:
        - harmonic: {index: [0,1], k: 100, req: 1.5}
        - ...
  - carbon_dioxide:
      structure:
        - O: [-1.162,0,0]
        - C: [0,0,0]
        - O: [1.162,0,0]
      bondlist:
        - harmonic: {index: [1,0], k: 8443, req: 1.162}
        - harmonic: {index: [1,2], k: 8443, req: 1.162}
        - harmonic_torsion: {index: [0,1,2], k: 451.9, aeq: 180}
      excluded_neighbours: 2 # generates an exclusionlist as shown below
      exclusionlist: [ [0,1], [1,2], [0,2] ] # redundant in this topology
  - ...

Structure Loading Policies

When giving structures using the structure keyword, the following policies apply:

  • structure can be a file name: file.@ where @=xyz|pqr|aam|gro
  • structure can be an array of atom names and their positions: - Mg: [2.0,0.1,2.0]
  • structure can be a FASTA sequence: {fasta: AAAAAAAK, k: 2.0; req: 7.0} which generates a linear chain of harmonically connected atoms. FASTA letters are translated into three letter residue names which must be defined in atomlist. Special letters: n=NTR, c=CTR, a=ANK. Instead of a sequence, fasta may be a filename from which the first sequence is extracted. The filename must end with .fasta.
  • Radii in files are ignored; atomlist definitions are used. A warning is issued if radii/charges differ in files and atomlist.
  • By default, charges in files are used; atomlist definitions are ignored. Use keepcharges=False to override.
  • If the structure file contains box size information, this will be ignored.

Nonbonded Interaction Exclusion

Some nonbonded interactions between atoms within a molecule may be excluded in the topology. Force fields almost always exclude nonbonded interactions between directly bonded atoms. However other nonbonded interactions may be excluded as well; refer to your force field parametrization. If a molecule contains overlapping hard spheres, e.g., if the bond length is shorter than the spheresʼ diameter, it is necessary to exclude corresponding nonbonded interactions to avoid infinite energies.

The excluded nonbonded interactions can be given as an explicit list of atom pairs excludelist, or they can be deduced from the moleculeʼs topology using the excluded_neighbours=n option: If the atoms are n or less bonds apart from each other in the molecule, the nonbonded interactions between them are excluded. Both options excluded_neighbours and exclusionlist can be used together making a union.

Initial Configuration

Upon starting a simulation, an initial configuration is required and must be specified in the section insertmolecules as a list of valid molecule names. Molecules are inserted in the given order and may be inactive, meaning that they are not present in the simulation cell, but available as a reservoir for e.g. grand canonical moves. If a group is marked atomic, its atoms are inserted N times.

Example:

insertmolecules:
  - salt:  { molarity: 0.1 }
  - water: { N: 256, inactive: 2 }

The following keywords for each molecule type are available:

insertmolecules Description
N Number of molecules to insert
molarity Insert molecules to reach molarity
inactive Number of inserted molecules to deactivate; set to true for all
positions Load positions from file (aam, pqr, xyz)
translate=[0,0,0] Displace loaded positions with vector

A filename with positions for the N molecules can be given with positions. The file must contain exactly N-times molecular positions that must all fit within the simulation box. Only positions from the file are copied; all other information is ignored.

For implicit molecules, only N should be given and the molecules are never inserted into the simulation box.

The molarity keyword is an alternative to N and uses the initial volume to calculate the number of molecules to insert. N and molarity are mutually exclusive.

Overlap Check

Random insertion is repeated until there is no overlap with the simulation container boundaries. Overlap between particles is ignored and for i.e. hard-sphere potentials the initial energy may be infinite.

Equilibrium Reactions

Faunus supports density fluctuations, coupled to chemical equilibria with explicit and/or implicit particles via their chemical potentials as defined in the reactionlist detailed below, as well as in atomlist and moleculelist. The level of flexibility is high and reactions can be freely composed.

The move involves deletion and insertion of reactants and products and it is therefore important that simulations are started with a sufficiently high number of initial molecules in insertmolecules. If not, the rcmc move will attempt to issue warnings with suggestions how to fix it.

Reaction format

The initial key describes a transformation of reactants (left of =) into products (right of =) that may be a mix of atomic and molecular species.

  • all species, +, and = must be surrounded by white-space
  • atom and molecule names cannot overlap
  • species can be repeated to match the desired stoichiometry, e.g. A + A = C

Available keywords:

reactionlist Description
lnK/pK Molar equilibrium constant either as $\ln K$ or $-\log_{10}(K)$
neutral=false If true, only neutral molecules participate in the reaction

The neutral keyword is needed for molecular groups containing titratable atoms. If neutral is set to true, the activity of the neutral molecule should be specified in moleculelist.

Example: Grand Canonical Salt Particles

This illustrates how to maintain constant chemical potential of salt ions:

atomlist:
  - na: {q: 1.0, ...}  # note that atom names must differ
  - cl: {q: -1.0, ...} # from molecule names
moleculelist:
  - Na+: {atoms: [na], atomic: true, activity: 0.1}
  - Cl-: {atoms: [cl], atomic: true, activity: 0.1}
reactionlist:
  - "= Na+ + Cl+": {} # note: molecules, not atoms
moves:
  - rcmc: {} # activate speciation move

The same setup can be used also for molecular molecules, i.e. molecules with atomic: false.

Example: Acid-base titration with implicit protons

An implicit atomic reactant or product is included in the reaction but not explicitly in the simulation cell. Common use-cases are acid-base equilibria where the proton concentration is often very low:

atomlist:
  - H+: {implicit: true, activity: 0.00001} # pH 5
  - COO-: {q: -1.0, ...}
  - COOH: {q: 0.0, ...}
reactionlist:
  - "COOH = COO- + H+": {pK: 4.8} # not electroneutral!

where we set pK equal to the pKa: $$ K_a = \frac{ a_{\mathrm{COO^-}} a_{\mathrm{H^+}} }{ a_{\mathrm{COOH}} }. $$ To simulate at a given constant pH, H+ is specified as an implicit atom of activity $10^{-\mathrm{pH}}$ and the equilibrium is modified accordingly (in this case $K$ is divided by $a_{\mathrm{H^+}}$). It is important to note that this reaction violates electroneutrality and should be used only with Hamiltonians where this is allowed. This could for example be in systems with salt screened Yukawa interactions.

Example: Acid-base titration coupled with Grand Canonical Salt

To respect electroneutrality when swapping species, we can associate the titration move with an artificial insertion or deletion of salt ions. These ions should be present under constant chemical potential and we therefore couple to a grand canonical salt bath:

atomlist:
  - H+: {implicit: true, activity: 0.00001} # pH 5
  - COO-: {q: -1.0, ...}
  - COOH: {q: 0.0, ...}
  - na: {q: 1.0, ...}
  - cl: {q: -1.0, ...}
moleculelist:
  - Na+: {atoms: [na], atomic: true, activity: 0.1}
  - Cl-: {atoms: [cl], atomic: true, activity: 0.1}
reactionlist:
  - "COOH + Cl- = COO- + H+": {pK: 4.8} # electroneutral!
  - "COOH = Na+ + COO- + H+": {pK: 4.8} # electroneutral!
  - "= Na+ + Cl-": {} # grand canonical salt

For the first reaction, $K$ is divided by both $a_{\mathrm{H^+}}$ and $a_{\mathrm{Cl^-}}$, so that the final equilibrium constant used by the speciation move is $$ K’ = \frac{K_a}{a_{ \mathrm{H^+} } a_{ \mathrm{Cl^-} } } = \frac{ a_{\mathrm{COO^-}} }{ a_{\mathrm{COOH}} a_{\mathrm{Cl^-}} }. $$ In an ideal system, the involvement of Na or Cl in the acid-base reaction is inconsequential for the equilibrium, since the Grand Canonical ensemble ensures constant salt activity.

Example: Precipitation of Calcium Hydroxide using implicit molecules

Here we introduce a solid phase of Ca(OH)₂ and its solubility product to predict the amount of dissolved calcium and hydroxide ions. Note that we start from an empty simulation box (both ions are inactive) and the solid phase is treated implicitly, i.e. it never inters the simulation box. If a forward reaction is made, one implicit molecule (out 200 in total) will be converted into explicit molecules. Once all 200 have been consumed, only backward reactions are possible. Additional, coupled reactions can be introduced to study complex equilibrium systems under influence of intermolecular interactions.

moleculelist:
    - Ca(OH)2: {implicit: true} # this molecule is implicit
    - Ca++: {atoms: [ca++], atomic: true}
    - OH-: {atoms: [oh-], atomic: true}
insertmolecules:
    - Ca++: {N: 200, inactive: true}
    - OH-: {N: 400, inactive: true}
    - Ca(OH)2: {N: 200} # not actually inserted!
reactionlist:
    - "Ca(OH)2 = Ca++ + OH- + OH-": {pK: 5.19}

Example: Swapping between molecular conformations

The following can be used to alternate between different molecular conformations. When swapping, the mass center position and orientation are randomly generated.

moleculelist:
  - A: {structure: "conformationA.xyz", ...}
  - B: {structure: "conformationB.xyz", ...}
reactionlist:
  - "A = B": {lnK: 0.69} # K=2, "B" twice as likely as "A"

Internal degrees of freedom (experimental)

If a fluctuating molecule has internal degreees of freedom, the internal bond energy is included as a bias so that the internal state does not affect the acceptance. To disable this behavior, a minor code modification is currently required (see MolecularGroupDeActivator::apply_bond_bias).