MW System Configuration
=======================

The goal of this section is to detail the syntax of the configuration input file. In MW, the system is configured by the user via an ASCII input
file *runtime.inpt*. Keywords are used to define a specific parameter or to introduce a block of parameters. 

First some global remarks about the configuration file syntax:

* Keywords are all in lower cases.
* One or more space/tab characters indicate the end of a keyword/parameter. This means that string parameters cannot have space in them, use underscore instead.
* There must be at most one keyword per line.
* The parameters associated with the keyword must be on the same line as
  the keyword.
* Blank lines and indentation are allowed in the configuration file,
  but not mandatory. They should be used to increase the readability.
* Comments are introduced by a '#' or '!' character and continue until
  the end of the line.
* The order of the keywords does not matter, however keywords
  belonging to the same block must be grouped together.
* If not stated otherwise, the parameters are given in atomic units.
* Some parameters have default values indicated in the following with `# default`. If the keyword is not specified, default values are used. Optional parameters are indicated in square brackets.


[Top level parameter keywords](#toplevelparameter)
* [num_steps](#num_steps)
* [timestep](#timestep)
* [temperature](#temperature)
* [num_pbc](#num_pbc)

[Top level block keywords](#toplevelblock)
* [thermostat](#thermostat)
* [barostat](#barostat)
* [velocity](#velocity)
* [species](#species)
   - [species_type](#species_type)
   - [dipoles_minimization](#dipoles_minimization)
   - [radius_minimization](#radius_minimization)
* [molecules](#molecules)
   - [constraint](#constraint)
   - [harmonic_bond](#harmonic_bond)
   - [harmonic_angle](#harmonic_angle)
   - [dihedral](#dihedral)
   - [improper](#improper)
* [electrodes](#electrodes)
   - [electrode_type](#electrode_type)
   - [electrode_charges](#electrode_charges)
   - [charge_neutrality](#charge_neutrality)
   - [compute_force](#compute_force)
* [interactions](#interactions)
   - [coulomb](#coulomb)
   - [lennard-jones](#lennard-jones)
   - [fumi-tosi](#fumi-tosi)
   - [xft](#xft)
   - [daim](#daim)
   - [damping](#damping)
   - [steele](#steele)
* [external_field](#external_field)
* [plumed](#plumed)
* [output](#output)

<a name="toplevelparameter">

Top level parameter keywords
----------------------------
</a>

The top-level keywords are used to define global system parameters
which do not fit in any specific category.

<a name="num_steps">

### num_steps

</a>

```
num_steps *int*
```

The *num_steps* keyword defines the number of timesteps in the
dynamics.

<a name="timestep">

### timestep

</a>

```
timestep *real*
```

The *timestep* keyword defines the duration of a timestep when the
system evolves. The value of the timestep is given in atomic
units, 1 a.u. = 2.418884·10<sup>−17</sup> s =  0.02418884 fs. A typical value of 41.341 corresponds to 1 fs.

<a name="temperature">

### temperature

</a>

```
temperature  *real*
```
The *temperature* keyword defines the system temperature. Its
value is given in K. This value is used when velocities need to be
initialized, when temperature rescaling is requested or when a 
thermostat is used. It is usually different from the actual system
temperature at a given time step derived from the kinetic energy of the
mobile particles.
A warning is raised in the *run.out* when the individual temperature of a particle exceeds a certain threshold.

<a name="num_pbc">

### num_pbc

</a>

```
num_pbc  *int*
```
The *num_pbc* keyword defines the number of periodic boundary conditions in the system. The value can be either *2* (PBC in *x* and *y* directions) or *3* (PBC in *x*, *y* and *z* directions).

<a name="toplevelblock">

Top-level block keywords
------------------------
</a>

<a name="thermostat">

### thermostat

</a>

```
thermostat
   chain_length    *int*      # default 5
   relaxation_time *real*     # default 4134.1
   max_iteration   *int*      # default 50
   tolerance       *real*     # default 1.0e-15
```
In order to simulate a system in the NVT ensemble, Nosé–Hoover chain
thermostat is implemented as described [here](NVT-ensemble). The
*chain_length* parameter defines the number of thermostats in the
chain, the first thermostat is coupled to the system and the second
thermostat, the second thermostat is coupled to the first and third
thermostat, and so on.

The *relaxation_time* parameter is the characteristic period of the thermostat given in atomic unit (1 a.u. = 2.418884·10<sup>−17</sup> s). 

Since the thermostat and the system are coupled, an iterative procedure is used to rescale the velocities of the particles and of the thermostats. The *max_iteration* and *tolerance* parameters are used to control this iterative process. 
The iterative process stops when $`\sum{m_i|v_i^{k}-v_i^{k-1}|^2} < \text{tol} N_{dof} k_b T`$.

The presence of the *thermostat* block indicates that the user is running a simulation in the NVT ensemble. In order to perform an NVE simulation, this block must be removed (or commented out) entirely.

</a>

<a name="barostat">

### barostat

</a>

```
barostat
   pressure        *real*     # default 0.0
   chain_length    *int*      # default 5
   relaxation_time *real*     # default 20670.5
```
This keyword is used to simulate a system in the NPT ensemble (see [here](NPT-ensemble) for a description of the implementation. Note that it is compulsory to define a thermostat for the ionic degrees of freedom (keyword [thermostat](#thermostat)). The NPT ensemble can be used only for 3D PBCs and without any electrode (in such case an alternative is to use the *piston* parameter in [electrode_type](#electrode_type)).

The *pressure* parameter is the external pressure which is applied to the system, given in atomic unit (1 a.u. = 2.942191·10<sup>13</sup> Pa).

In the method, the barostat degree of freedom is also thermostated. The
*chain_length* parameter defines the number of thermostats in the
chain, the first thermostat is coupled to the barostat and the second
thermostat, the second thermostat is coupled to the first and third
thermostat, and so on.

The *relaxation_time* parameter is the characteristic period of the thermostat given in atomic unit (1 a.u. = 2.418884·10<sup>−17</sup> s). 

Since the barostat and the system are coupled, an iterative procedure is used to rescale the velocities of the particles and of the thermostats/barostat. The parameters used in this iterative procedure are specified using the [thermostat](#thermostat) keyword.

<a name="velocity">

### velocity

</a>

```
velocity
   create          *true/false*         # default false if velocities given in *data.inpt*, true otherwise
   com_threshold   *real*               # default 1.0e+6
   scale           *string* *int*       # default global 100
```
The velocity block is used to define the following actions:

* *create*: if velocities are not provided in the *data.inpt* file, then these are always sampled from a pseudo-random normal distribution with fixed seed and rescaled at the system temperature, regardless of the value of the keyword. If velocities are given in the *data.inpt* file, then the keyword *create* can be used to decide whether to resample the velocities at the beginning of the simulation from a pseudo-random normal distribution at the given temperature or not. 
* *com_threshold*: the parameter is a threshold value which triggers a reset of the centre of mass of all mobile particles if the kinetic energy of the centre of mass exceeds the threshold.
* *scale*: the *scale* parameters determine the method used to scale velocities and at which period (in number of steps) this scaling occurs, if any. The scale method keywords are:
    + *global*: all species are scaled together
    + *species*: species types are scaled independently from each other
    + *molecules_global*: all species are scaled together, species belonging to the same molecule are seen as one species type
    + *molecules_independent*: species types are scaled independently from each other, species belonging to the same molecule are seen as one species type.

If *scale* keyword is not given, no scaling occurs except at the beginning of the simulation if velocities are created. In this case the method used is *global*. If *scale* is given without further specifications, default method is *global* and the scaling frequency is *100*.

<a name="species">

### species

</a>

Low-level blocks:
* [species_type](#species_type)
* [dipoles_minimization](#dipoles_minimization)
* [radius_minimization](#radius_minimization)

```
species
   species_type
       name   *string*  
       count  *int*                       # default 0
       mass   *real*                      # default 0.0
       charge *string* [*real*] [*real*]  # default neutral
       polarizability *real*              # default 0.0
       deformability *real* *real* *real* # default 0.0 0.0 0.0
       mobile *true/false*                # default true

    repeat species_type block for each species
    ...

   dipoles_minimization *string* [*string*] [*real*] [*int*]   # default cg 1.0e-12 100 
   radius_minimization *string* [*real*] [*int*]   # default cg 1.0e-12 100 
```

The *species* block is used to define species which are simulated in the system. 

<a name="species_type">

#### species_type

</a>

For each species, one must define a *species_type* block.

The *name* parameter must be a unique identifier for this species. It can be made of up to 8 ASCII characters but cannot contain white spaces, '#' or '!'.

The *count* parameter is the number of species of this type in the system.

The *mass* parameter is the mass of this species type in the unified atomic mass unit (u).

The *charge* parameter sets the charge type and its parameters for this species type. It can be:

* `charge neutral`: defines a species with no electrostatic interaction
* `charge point *real*`: defines a species with a point-charge distribution. The *real* parameter is the magnitude of the charge in the atomic unit (e)
* `charge gaussian *real* *real*`: defines a species with a gaussian-charge distribution. The first *real* parameter is the $`\eta`$ width parameter, the second *real* parameter is the magnitude of the charge in the atomic unit (e).

The *polarizability* parameter sets the polarizability of the species as described in the [force field](force-field) section.

Three *deformability* parameters defines the deformability of the species as described in the [force field](force-field) section. 

It is not possible to specify *polarizability* and *deformability* when there are electrodes in the simulation.
The *mobile* parameter determines if particles of this species type are mobile during the dynamics evolution.

The species must be given in the same order as they are in the data file *data.inpt*. Species that do not belong to electrodes must come first. They must have a point charge distribution or be of a neutral type.

Species in the electrodes must be given last. They must have a Gaussian charge distribution. They cannot be mobile.

<a name="dipoles_minimization">

#### dipoles_minimization

</a>

If at least one atomic species is polarizable, an algorithm for dipoles computation must be specified.

Two methods are available to compute the ionic dipoles $`\bm{\mu}`$ at each time step. 
The first solves directly the equation $`\bm{\mu} = \alpha\bm{E}`$ computing the electric field on one ions due to the others ($`\alpha`$ is the polarizability of the ion). This problem can be solved *via* conjugate gradient minimization using
```
dipoles_minimization cg *real* *int*          # default 1.0e-12 100
```
where the first parameter is the tolerance criterion for the convergence and the second parameter is the maximum number of iterations allowed. 
The iterative process stops when $`|\text{residual}| < \text{tol} * \sqrt{3N_{ions}}`$ where $`N_{ions}`$ refers only to the polarizable ones.
The same equation can be solved by direct matrix inversion using 
```
dipoles_minimization matrix_inversion
```
However, this option is not recommended as the matrix has to be computed and inverted at each timestep, which is usually not efficient.

The second method is the mass-zero constrained molecular dynamics (MaZe), denoted by the keyword *maze*. The dipoles are computed solving the mass-zero constrained equations obtained by the adiabatic condition $`\bm{\nabla}_{\bm{\mu}}U = 0`$ using algorithms based on SHAKE.
For MaZe algorithm using an iterative SHAKE algorithm use
```
dipoles_minimization maze iterative_shake *real* *int* [*int*]   # default 1.0e-12 100 0
```
where the first parameter is the tolerance criterion for the convergence and the second parameter is the maximum number of iterations allowed. 
The third parameter is called *nblocks* and is optional. It controls the level of approximation for the SHAKE matrix. 
If 0, not given or equal to the number of electrode atoms, the matrix is approximated by using only its diagonal elements. 
If greater than 0, the SHAKE matrix is approximated by using *nblocks* block diagonal elements. It should be chosen such that the number of electrode atoms is a multiple of *nblocks*.
The iterative process stops when $`\text{max}(|\bm{\nabla}_{\bm{\mu}} U|) < \text{tol}`$.

The shake equation can be solved also by direct inversion of the shake matrix using
```
dipoles_minimization maze inversion 
```
However, this last option is not recommended as the matrix has to be computed and inverted at each timestep, which is usually not efficient.

It is not possible to specify a method for dipoles computation when electrodes are present in the simulation.

<a name="radius_minimization">

#### radius_minimization

</a>

If at least one atomic species is deformable (as in the Dipole Aspherical Ion Model, DAIM), an algorithm for radii computation must be specified.

The only algorithm available is the nonlinear conjugate gradient which can be specified using
```
radius_minimization cg *real* *int*          # default 1.0e-12 100
```
where the first parameter is the tolerance criterion for the convergence and the second parameter is the maximum number of iterations allowed. 
The iterative process stops when $`|\text{residual}| < \text{tol} * \sqrt{3N_{ions}}`$ where $`N_{ions}`$ refers only to the deformable ones.

<a name="molecules">

### molecules

</a>

```
molecules
   molecule_type
     name *string*
     count *int*                                  # default 0
     sites *string* *string* ...    
   
     constraints_algorithm *string* *parameters*  # default rattle
     constraint *string* *string* *real*
     constraint *string* *string* *real*
     ... 

     harmonic_bond *string* *string* *real* *real*
     harmonic_bond *string* *string* *real* *real*
     ...

     harmonic_angle *string* *string* *string* *real* *real*
     harmonic_angle *string* *string* *string* *real* *real*
     ...

     dihedral *string* *string* *string* *string* *real* *real* *real* *real*
     dihedral *string* *string* *string* *string* *real* *real* *real* *real*
     ...

     improper *string* *string* *string* *string* *real* *real* *real* *real*
     improper *string* *string* *string* *string* *real* *real* *real* *real*
     ...

   repeat molecule_type block for each molecule
   ...
```

The *molecules* block introduces the definition of one or more molecule types. Each molecule type parameters are set in a *molecule_type* block.

The *name* parameter must be a unique identifier for this molecule type. It can be made of up to 8 ASCII characters but cannot contain white spaces, '#' or '!'.

The *count* parameter is the number of molecules of this type. It should be the same value as the *count* parameter for each site species.

The *sites* parameter defines the species which compose the molecule. Each species name must correspond to a unique species type defined in the *species* block. Because runtime checks and mappings are performed during the configuration parsing, the *species* block must be defined before the *molecules* block. There can be at most 59 sites per molecule.

<a name="constraint">

If the molecule has rigid constraints, then the parameters *constraints_algorithm* and *constraint* must be used to define them. The *constraints_algorithm* defines the algorithm used to solve the constraints. Currently, two algorithms are supported: *rattle* and *linear*. *rattle* is used to resolve generic rigid bond constraints while *linear* is used for 3 sites linear molecules. The two forms of *constraints_algorithm* parameters are
```
constraints_algorithm rattle *real* *int*   # default 1.0e-6 100
```
for the *rattle* algorithm where the first parameter corresponds to the tolerance on constraint realisation and the second parameter is the maximum number of iterations.  For more details on the meaning of these two parameters please refer to the publications listed [here](time-integration).

```
constraints_algorithm linear *string* *string*
```
for the *linear* algorithm where the two parameters correspond to the two sites used as the basis of the linear molecule. The algorithm implementation is described in the pdf file *./doc/theory/RATTLE-linear-triatomic-molecules/rattle_linear_triatomic_molecule.pdf*.

The *constraint* keyword is used to set the bond length between two sites of the molecule:
```
constraint *string* *string* *real*
```
The first and second parameters correspond the sites names involved in this link and the third parameter is the rigid bond length. All constraints must be defined consecutively. Constraints must not be redundant. For linear molecules, exactly 3 constraints must be defined.

Intra-molecular potentials can be defined with the *harmonic_bond*, *harmonic_angle*, *dihedral* and *improper* keywords.

<a name="harmonic_bond">

The *harmonic_bond* keyword defines a harmonic potential between two sites of the molecule of the form 
```math
V(r_{ij}) = k_0(|r_{ij}| - r_0)^2
``` 
with the following syntax:
```
harmonic_bond *string* *string* *real* *real*
```
where the first and second parameters correspond to the identifier of the sites involved in the harmonic bond. The third parameter is the strength of the harmonic bond (in atomic units: $`E_h / a_0`$) and the last parameter is the length. All harmonic bonds must be defined consecutively.

<a name="harmonic_angle">

The *harmonic_angle* keyword defines a harmonic potential between three sites of the molecule of the form
```math
V(r_i,r_j,r_k) = k_0(\theta_{ijk} - \theta_0)^2
```
with the following syntax:
```
harmonic_angle *string* *string* *string* *real* *real*
```
where the first, second and third parameters correspond to the identifier of the sites involved in the harmonic angle. The sites of the $`\hat{ijk}`$ angle should be given in the order $`j i k`$ (the central site first). The fourth parameter is the strength of the harmonic angle (in atomic units: $`E_h / rad`$) and the last parameter is the angle at equilibrium (in rad). All harmonic angles must be defined consecutively.

<a name="dihedral">

The *dihedral* keyword defines a torsion potential between four sites of the molecule of the form 
```math
V(r_i,r_j,r_k,r_l) = 0.5*(v_1(1+cos(\phi_{ijkl}))+v_2(1-cos(2\phi_{ijkl}))+v_3(1+cos(3\phi_{ijkl}))+v_4*(1-cos(4\phi_{ijkl})))
```
with the following syntax:
```
dihedral *string* *string* *string* *string* *real* *real* *real* *real*
```
where the first, second, third and fourth parameters correspond to the identifier of the sites involved in the dihedral angle. The sites of the $`ijkl`$ angle should be given in the order $`i j k l`$ (in consecutive order). The fifth to eighth parameters are the $`v_1`$, $`v_2`$, $`v_3`$, $`v_4`$ parameters (in atomic units: $`E_h`$). All dihedrals must be defined consecutively.

<a name="improper">

The *improper* keyword defines an improper dihedral potential between four sites of the molecule of the form 
```math
V(r_i,r_j,r_k,r_l) = 0.5*(v_1(1+cos(\phi_{ijkl}))+v_2(1-cos(2\phi_{ijkl}))+v_3(1+cos(3\phi_{ijkl}))+v_4*(1-cos(4\phi_{ijkl})))
```
with the following syntax:
```
improper *string* *string* *string* *string* *real* *real* *real* *real*
```
where the first, second, third and fourth parameters correspond to the identifier of the sites involved in the improper angle. The sites of the $`ijkl`$ angle should be given with the central atom third. The fifth to eighth parameters are the $`v_1`$, $`v_2`$, $`v_3`$, $`v_4`$ parameters (in atomic units: $`E_h`$). All impropers must be defined consecutively.

All intramolecular interactions are excluded for 1-2 and 1-3 bonded sites. The 1-4 intramolecular interactions are scaled by a factor of 0.5.

<a name="electrodes">

### electrodes

</a>

Low-level blocks:
* [electrode_type](#electrode_type)
* [electrode_charges](#electrode_charges)
* [charge_neutrality](#charge_neutrality)
* [compute_force](#compute_force)

```
electrodes

   electrode_type
      name *string*
      species *string*                           
      potential *real*                           # default 0.0
      piston *real* *int*                        # default 0.0 1
      thomas_fermi_length *real*                 # default 0.0
      voronoi_volume *real*                             

   repeat electrode_type block for each electrode
   ...

   electrode_charges *string* [*string*] [*real*] [*int*]   # default cg 1.0e-12 100 
   charge_neutrality *true/false* [*global/electrodes*]     # default false global
   compute_force *true/false*       # default false
```

The *electrodes* block is used to define electrodes to which species with a Gaussian charge distribution for electrostatic interactions are assigned.
They cannot be defined when polarizable species are present in the simulation.

<a name="electrode_type">

#### electrode_type

</a>

For each electrode, the parameters are defined in a *electrode_type* block.

```
electrode_type
   name      *string*
   species   *species*
   potential *real*                           # default 0.0
   piston *real* *int*                        # default 0.0 1
   thomas_fermi_length *real*                 # default 0.0
   voronoi_volume *real*
```

The *name* parameter must be a unique identifier for this electrode species. It can be made of up to 8 ASCII characters but cannot contain white spaces, '#' or '!'.

The *species* parameter is the name of the species used to fill this electrode. *WARNING* Each electrode can only be formed of exactly one species type. To simulate multiple species electrode (for example a magnesium oxide electrode MgO, one would have to define two electrodes per MgO electrode. To move them with a piston, it is *not* possible using directly MW, but one can use the [python interface to MW](python-interface) to do so, see the example `./example/python-piston-mgo`)

The *potential* parameter sets the constant potential for this electrode type.

The *piston* parameter triggers an NPT-like simulation by using the electrodes as rigid pistons, as described [here](NPzT-ensemble). The electrodes are evolved by adding a term equal to $`P \times S`$ where $`S`$ is the electrode section and $`P`$ is the target pressure. The first parameter is the target pressure. The second parameter is an integer that gives the direction of the applied force (should be $`\pm 1`$). This option automatically sets *compute_force* to *true*.

The *thomas_fermi_length* parameter sets the Thomas-Fermi length for this electrode.

The *voronoi_volume* parameter sets the Voronoi volume for this electrode.

<a name="electrode_charges">

#### electrode_charges

</a>

The charges can either be kept constant with *constant_charge* keyword or evolved to keep the potential constant. 
For constant charges electrodes use
```
electrode_charges constant_charge
```
The charge used to initialize the electrode charges is the one of the associated species or the one given in the data file if present.

If the simulation is at constant potential, the *electrode_charges* keyword gives the algorithm used to evolve the charges on the electrodes. 
Two methods are available to compute the charges on electrodes at each time step. 
The first solves directly the equation $`\frac{\partial U}{\partial Q} = \Psi`$ using the electrode-electrode interaction matrix. This problem can be solved 
by inverting this matrix, using the keyword *matrix_inversion*, or by using a conjugate gradient with the keyword *cg*.
For constant potential electrodes with matrix inversion use
```
electrode_charges matrix_inversion
```
The first time a system is run, the matrix is computed and stored in the file *hessian_matrix.inpt*. If the file already exists, instead of being computed, the matrix is read from that file. A consistency check on the first atom position is made. 

For constant potential electrodes with conjugate gradient use
```
electrode_charges cg *real* *int*          # default 1.0e-12 100
```
where the first parameter is the tolerance criterion for the convergence and the second parameter is the maximum number of iterations allowed. 
The iterative process stops when $`|\text{residual}| < \text{tol} * \sqrt{N_{atoms}}`$.

The second method is the mass-zero constrained molecular dynamics (MaZe), denoted by the keyword *maze*. Several implementations are possible that solve the constraints
by inverting a matrix (*inversion*), using an iterative SHAKE procedure (*iterative_shake*) or by using a conjugate gradient procedure (*iterative_cg*).
For MaZe algorithm using the SHAKE matrix inversion use
```
electrode_charges maze inversion 
```
For MaZe algorithm using an iterative SHAKE algorithm use
```
electrode_charges maze iterative_shake *real* *int* [*int*]   # default 1.0e-12 100 0
```
where the first parameter is the tolerance criterion for the convergence and the second parameter is the maximum number of iterations allowed. 
The third parameter is called *nblocks* and is optional. It controls the level of approximation for the SHAKE matrix. 
If 0, not given or equal to the number of electrode atoms, the matrix is approximated by using only its diagonal elements. 
If greater than 0, the SHAKE matrix is approximated by using *nblocks* block diagonal elements. It should be chosen such that the number of electrode atoms is a multiple of *nblocks*.
The iterative process stops when $`\text{max}(|\frac{\partial U}{\partial Q} - \Psi|) < \text{tol}`$.
For MaZe algorithm using an iterative conjugate gradient algorithm use
```
electrode_charges maze iterative_cg *real* *int*   # default 1.0e-12 100
```
where the first parameter is the tolerance criterion for the convergence and the second parameter is the maximum number of iterations allowed. 

The first time a system is run, the matrix is computed and stored in the file *maze_matrix.inpt*. If the file already exists, instead of being computed, the matrix is read from that file. A consistency check on the first atom position is made. 

<a name="charge_neutrality">

#### charge_neutrality

</a>

The *charge_neutrality* option si used to perform simulation at constant potential under the so-called electroneutrality constraint. If the second parameter is set to *global* (default), the total charge on the electrodes is constrained to be the opposite of the electrolyte charge, so that the system is globally neutral. If the second parameter is set to *electrodes*, the total charge is constrained to be zero only for the electrode system. All algorithms for solving the electrode charges can be used (constrained conjugate gradient, MaZe with additional conditions -- only the *inversion* case -- or modified matrix inversion). If the algorithm chosen is `maze` the parameter $`\kappa`$ (see the page [Constant Potential Simulation](constant-potential) and references therein for more details) can also be specified after the second keyword (*global* or *electrodes*). The syntax for this option is then
```
charge_neutrality true/false [global/electrodes [*real*]]   # default false  
```

<a name="compute_force">

#### compute_force

</a>

The *compute_force* option allows to compute the total force on the electrodes.

<a name="interactions">

### interactions

</a>

Interatomic interactions are defined in the *interactions* block. Low-level blocks are:
- [coulomb](#coulomb)
- [lennard-jones](#lennard-jones)
- [fumi-tosi](#fumi-tosi)
- [xft](#xft)
- [daim](#daim)
- [damping](#damping)
- [steele](#steele)

```
interactions
   coulomb
      coulomb_rcut *real*     # default 0.0
      coulomb_rtol *real*     # default 1.0e-15
      coulomb_ktol *real*     # default 1.0e-15
   lennard-jones 
      lj_rcut *real*
      [lj_3D_tail_correction]
      lj_pair *string* *string* *real* *real*
      lj_pair *string* *string* *real* *real*
      ...
      [lj_rule *string*]
   fumi-tosi
      ft_rcut *real*
      [ft_3D_pressure_tail_correction]
      ft_pair *string* *string* *real* *real* *real* *real* *real* *real*
      ft_pair *string* *string* *real* *real* *real* *real* *real* *real*
      ...
      [ft_rule *string*]
   xft
      xft_rcut *real*
      [xft_3D_pressure_tail_correction]
      xft_pair *string* *string* *real* *real* *int* *real* *real* *real* *real* *real* *real*
      xft_pair *string* *string* *real* *real* *int* *real* *real* *real* *real* *real* *real*
      ...
      [xft_rule *string*]
   daim
      daim_rcut *real*
      [daim_3D_pressure_tail_correction]
      daim_pair *string* *string* *real* *real* *real* *real* *real* *real* *real* *real* *real* *real*
      daim_pair *string* *string* *real* *real* *real* *real* *real* *real* *real* *real* *real* *real*
      ...
   damping
      tt_pair *string* *string* *real* *int* *real*
      tt_pair *string* *string* *real* *int* *real*
      ... 
   steele *int*
      steele_rcut *real*
      wall *string* *string* *real* *real* *real* *real*
      wall *string* *string* *real* *real* *real* *real*
      ...
```

<a name="coulomb">

#### coulomb

</a>

Electrostatic interactions are treated with the Ewald summation formalism to properly handle the long-range interactions in periodic boundary conditions. Ewald parameters are defined in the *coulomb* block. 

```
coulomb
   coulomb_rcut *real*     # default 0.0
   coulomb_rtol *real*     # default 1.0e-15
   coulomb_ktol *real*     # default 1.0e-15
```

*coulomb_rcut* sets the cut-off distance for real space interactions. Since the implementation uses the minimum image distance convention this value should be less than half the minimum dimension of the box. An error is raised otherwise.

*coulomb_rtol* is the magnitude below which real-space term is not included in the summation. It is used in combination with *coulomb_rcut* to compute the $`\alpha`$ which splits computation between real space and reciprocal space.

*coulomb_ktol*  is the magnitude below which reciprocal-space term is not included in the summation. It is used to compute the number of k-points for the reciprocal space sum.

<a name="lennard-jones">

#### lennard-jones

</a>

Lennard-Jones interaction potential can be used to model the interatomic interactions. This potential is defined as 
```math
V(r_{ij}) = 4 \epsilon_{ij} \left( \frac{\sigma_{ij}^{12}}{r_{ij}^{12}} - \frac{\sigma_{ij}^6}{r_{ij}^{6}} \right)
```
The parameters of these interactions are defined in the *lennard-jones* block.

```
lennard-jones
   lj_rcut *real*
   [lj_3D_tail_correction]
   lj_pair *string* *string* *real* *real*
   lj_pair *string* *string* *real* *real*
   ...
   [lj_rule *string*]
```

*lj_rcut* sets the cut-off distance (in atomic units $`a_0`$) for all Lennard-Jones interactions. Since the implementation uses the minimum image distance convention this value should be less than half the minimum dimension of the box. An error is raised otherwise.

By default the Lennard-Jones potential implemented is shifted so that its value at the cut-off is zero. If the *lj_3D_tail_correction* optional keyword is used, the potential is not shifted anymore and a long-range tail correction is added to the energy and the virial. Since this correction assumes the pair correlation function $`g(r)`$ is unity after the cutoff, this should be used only in a fluid 3D system with no interfaces.
These modifications have no impact on the forces applied on the system.

For each pair of species interacting via a Lennard-Jones potential, one should provide *lj_pair* parameters. The first and second parameters correspond to the name of the species involved in the interaction. The names must belong to one the *species_type* definitions. The third parameter is the $`\epsilon`$ value in kJ/mol and the fourth parameter is $`\sigma`$ in angstrom. Be aware that these parameters are **not** in atomic units. If some pairs are missing, they do not interact with this potential by default.

If given, the optional keyword *lj_rule* calculates cross Lennard-Jones parameters ($`i-j`$) from direct parameters ($`i-i`$ and $`j-j`$) using mixing rules. The keyword should be placed after the necessary direct parameters have been defined with *lj_pair*. Any following *lj_pair* keyword will overwrite the calculated cross-term.
The mixing rule is chosen with the first parameter among the following strings: "Lorentz_Berthelot". The Lorentz-Berthelot rule states that 
```math
\epsilon_{ij} = \sqrt{\epsilon_{ii}\epsilon_{jj}} \\
\sigma_{ij} = (\sigma_{ii} + \sigma_{jj})/2
```

<a name="fumi-tosi">

#### fumi-tosi

</a>

The Fumi-Tosi interaction potential can be used to model the interatomic interactions. This potential is defined as 
```math
V_{ij}(r_{ij}) = B_{ij} \exp(-\eta _{ij} r_{ij}) - f_{dd}^{ij}(r_{ij}) \frac{C_{ij}}{r_{ij}^{6}} - f_{dq}^{ij}(r_{ij}) \frac{D_{ij}}{r_{ij}^{8}}
```
where
```math
f_{dd}^{ij}(r_{ij}) = 1 - \exp(-d_{dd}^{ij} r_{ij}) \sum_{k=0}^6 \frac{(d_{dd}^{ij} r_{ij})^k}{k!}
```
and 
```math
f_{dq}^{ij}(r_{ij}) = 1 - \exp(-d_{dq}^{ij} r_{ij}) \sum_{k=0}^8 \frac{(d_{dq}^{ij} r_{ij})^k}{k!}
```
The parameters of these interactions are defined in *fumi-tosi* block.

```
fumi-tosi
   ft_rcut *real*
   [ft_3D_pressure_tail_correction]
   ft_pair *string* *string* *real* *real* *real* *real* *real* *real*
   ft_pair *string* *string* *real* *real* *real* *real* *real* *real*
   ...
   [ft_rule *string*]
```

*ft_rcut* sets the cut-off distance for all Fumi-Tosi interactions. Since the implementation uses the minimum image distance convention this value should be less than half the minimum dimension of the box. An error is raised otherwise.

By default the Fumi-Tosi potential implemented is truncated at the cut-off distance. If the *ft_3D_pressure_tail_correction* optional keyword is used, a long-range tail correction is added to the virial for the dispersion terms. Since this correction assumes the pair correlation function $`g(r)`$ is unity after the cutoff, this should be used only in a fluid 3D system with no interfaces.
These modifications have no impact on the forces applied on the system.

For each pair of species interacting via a Fumi-Tosi potential, one should provide *ft_pair* parameters. The first and second parameters correspond to the name of the species involved in the interaction. The names must belong to one the *species_type* definitions. The third to eight parameters are respectively $`\eta_{ij}`$, $`B_{ij}`$, $`C_{ij}`$, $`D_{ij}`$, $`d_{dd}^{ij}`$, $`d_{dq}^{ij}`$ in atomic units.

If given, the optional keyword *ft_rule* calculates cross Fumi-Tosi parameters ($`i-j`$) from direct parameters ($`i-i`$ and $`j-j`$) using mixing rules. The keyword should be placed after the necessary direct parameters have been defined with *ft_pair*. Any following *ft_pair* keyword will overwrite the calculated cross-term.
The mixing rule is chosen with the first parameter among the following strings: "Mixing". The mixing rule is detailed in the [force field section](force-field) of the Wiki.

<a name="xft">

#### xft

</a>

The XFT interaction potential can be used to model the interatomic interactions. This potential is defined as 
```math
V_{ij}(r_{ij}) =\frac{B_{ij} \exp(-\eta_{ij}r_{ij})}{r_{ij}^{n}}  + B'_{ij} \exp(-\eta'_{ij} r^2_{ij}) - f_{dd}^{ij}(r_{ij}) \frac{C_{ij}}{r_{ij}^{6}} - f_{dq}^{ij}(r_{ij}) \frac{D_{ij}}{r_{ij}^{8}}
```
where $`f_{dd}^{ij}(r_{ij})`$ and $`f_{dq}^{ij}(r_{ij})`$ have the same definition as in the Fumi-Tosi potential.

The parameters of these interactions are defined in *xft* block.

```
xft
   xft_rcut *real*
   [xft_3D_pressure_tail_correction]
   xft_pair *string* *string* *real* *real* *int* *real* *real* *real* *real* *real* *real*
   xft_pair *string* *string* *real* *real* *int* *real* *real* *real* *real* *real* *real*
   ...
   [xft_rule *string*]
```

*xft_rcut* sets the cut-off distance for all XFT interactions. Since the implementation uses the minimum image distance convention this value should be less than half the minimum dimension of the box. An error is raised otherwise.

By default the XFT potential implemented is truncated at the cut-off distance. If the *xft_3D_pressure_tail_correction* optional keyword is used, a long-range tail correction is added to the virial for the dispersion terms. Since this correction assumes the pair correlation function $`g(r)`$ is unity after the cutoff, this should be used only in a fluid 3D system with no interfaces.
These modifications have no impact on the forces applied on the system.

For each pair of species interacting via a XFT potential, one should provide *xft_pair* parameters. The first and second parameters correspond to the name of the species involved in the interaction. The names must belong to one the *species_type* definitions. The third to eleven parameters are respectively $`\eta_{ij}`$, $`B_{ij}`$, $`n`$, $`\eta'_{ij}`$, $`B'_{ij}`$, $`C_{ij}`$, $`D_{ij}`$, $`d_{dd}^{ij}`$, $`d_{dq}^{ij}`$ in atomic units.

If given, the optional keyword *xft_rule* calculates cross XFT parameters ($`i-j`$) from direct parameters ($`i-i`$ and $`j-j`$) using mixing rules. The keyword should be placed after the necessary direct parameters have been defined with *xft_pair*. Any following *xft_pair* keyword will overwrite the calculated cross-term.
The mixing rule is chosen with the first parameter among the following strings: "Mixing". The mixing rule is detailed in the [force field section](force-field) of the Wiki.

<a name="daim">

#### daim

</a>

AIM potential can be used to model the interatomic interactions. This potential is defined as 
```math
V(r_{ij}) = B^1_{ij}\exp(-\eta^1_{ij} r_{ij}) + B^2_{ij}\exp(-\eta^2_{ij} \rho_{ij}) + B^3_{ij}\exp(-\eta^3_{ij} \rho_{ij}) - f_{dd}^{ij}(r_{ij})\frac{C_{ij}}{r_{ij}^6} - f_{dq}^{ij}(r_{ij})\frac{D_{ij}}{r_{ij}^8}
```
where $`f_{dd}^{ij}(r_{ij})`$ and $`f_{dq}^{ij}(r_{ij})`$ have the same definition as in the Fumi-Tosi potential.

The parameters of these interactions are defined in the *daim* block.

```
daim
   daim_rcut *real*
   [daim_3D_pressure_tail_correction]
   daim_pair *string* *string* *real* *real* *real* *real* *real* *real* *real* *real* *real* *real*
   daim_pair *string* *string* *real* *real* *real* *real* *real* *real* *real* *real* *real* *real*
   ...
```

*daim_rcut* sets the cut-off distance (in atomic units $`a_0`$) for all DAIM interactions. Since the implementation uses the minimum image distance convention this value should be less than half the minimum dimension of the box. An error is raised otherwise.

By default the DAIM potential implemented is truncated at the cut-off distance. If the *daim_3D_pressure_tail_correction* optional keyword is used, a long-range tail correction is added to the virial for the dispersion terms. Since this correction assumes the pair correlation function $`g(r)`$ is unity after the cutoff, this should be used only in a fluid 3D system with no interfaces.
These modifications have no impact on the forces applied on the system.

For each pair of species interacting via a DAIM potential, one should provide *daim_pair* parameters. The first and second parameters correspond to the name of the species involved in the interaction. The names must belong to one the *species_type* definitions. The third to thirteenth parameters are respectively $`\eta^1_{ij}`$, $`\eta^2_{ij}`$, $`\eta^3_{ij}`$, $`B^1_{ij}`$, $`B^2_{ij}`$, $`B^3_{ij}`$, $`C_{ij}`$, $`D_{ij}`$, $`d_{dd}^{ij}`$, $`d_{dq}^{ij}`$.

<a name="damping">

#### damping

</a>

This section provides the parameters of the Tang-Toennies damping functions for the charge-dipole interactions. These functions are defined as 
```math
f_n(r_{ij}) = 1 - c_n * exp(-b_n r_{ij}) \sum_{k=0}^{n} \frac{(b_n r_{ij})^k}{k!}
```
The parameters of these functions are defined in the *damping* block.

```
damping
   tt_pair *string* *string* *real* *int* *real*
   tt_pair *string* *string* *real* *int* *real*
   ...
```
For each pair of species, one should provide *tt_pair* parameters. The first and second parameters correspond to the name of the species involved. The first string corresponds to the atom with the charge whereas the second string corresponds to the atom with the dipole. As such, if two species bear both a charge and a dipole, the Tang-Toennies should be defined twice (once for the charge1-dipole2 and once for the charge2-dipole1). The names must belong to one the *species_type* definitions. The third to fifth parameters are respectively $`b_n`$, $`n`$, $`c_n`$.

<a name="steele">

#### steele

</a>

The Steele potential is a one-dimensional potential used to model non-structured walls. The potential is defined as
```math
V(|z_{i}-z_0|) = 2 \pi \rho_w \epsilon _{fw} \sigma _{fw} \Delta [2/5 \left(\frac{\sigma _{fw}}{z}\right) ^{10} - \left(\frac{\sigma _{fw}}{z}\right)^4 - \frac{\sigma_{fw}^4}{3 \Delta (z + 0.61 \Delta)^3}]
```
where $`z_i`$ is the distance to the wall. The wall should be represented by an *electrode_type* corresponding to a *species_type* containing only one atom and the coordinate of the atom will be taken as wall position $`z_0`$.

The parameters are defined in the *steele* block.

```
steele *int*
   steele_rcut *real*
   wall *string* *string* *real* *real* *real* *real*
   wall *string* *string* *real* *real* *real* *real*
   ...
```

The definition of the block *steele* requires an integer corresponding to the number of different walls (typically 1 or 2).
 
*steele_rcut* sets the cut-off distance (in atomic units $`a_0`$) for all Steele interactions. 

For each ion type interacting with the Steele wall, one should provide *wall* parameters. The first parameter is the electrode name corresponding to the wall, the second parameter is the species name interacting with the wall. The third to sixth parameter are respectively $`\rho_w`$ in angstroms$`^{-3}`$, $`\epsilon_{fw}`$ in kJ/mol, $`\sigma_{fw}`$ in angstroms and $`\Delta`$ in angstroms.

<a name="external_field">

### external_field

</a>

```
external_field
   type            *string*              
   direction       *real* *real* *real*   # default 0.0 0.0 0.0
   amplitude       *real*                 # default 0.0
```
The *external_field* block allows to use the finite field methods, either finite electric field or electric displacement:

* *type* sets the type of the applied field: *E* for electric field or *D* for electric displacement.
* *direction* sets the direction of the applied field with respect to the simulation box. It does not have to be a unit vector. The field will be applied only on the directions whose component is non zero.
* *amplitude* sets the quantity by which each component of the *direction* vector is multiplied to obtain the field.

Note that for the electric displacement, having a non zero *direction* and a 0.0 *amplitude* will apply a D=0 in the non zero components of *direction*, and this is not equivalent to not applying a field.

<a name="plumed">

### plumed

</a>

```
plumed *string*
```
MW can be run coupled with [Plumed](http://tcb.ucas.ac.cn/plumed2/user-doc/html/index.html) software to bias forces and many other things. MW has to be compiled with Plumed to enable this functionality (cf [installation instructions](install)). The first parameter is the name of the file that contains the Plumed instructions.

<a name="output">

### output

</a>

The *output* block defines which output files will be generated and the frequency at which they will be updated.

The *default* keyword is used to assign a default value for all output frequencies. Specific values can be set by using the corresponding keyword.

A value of 0 turns off output.

```
output
  default              *int*
  step                 *int*
  restart              *int* [alternate] 
  trajectories         *int*
  xyz                  *int*
  pdb                  *int*
  lammps               *int*
  charges              *int*
  total_charges        *int*
  energies             *int* [*string* *string* ...]
  forces               *int*
  elec_forces          *int*
  dipoles              *int*
  radius               *int*
  polarization         *int*
  potential_shift      *int*
  density              *num_bins* *frequency* [*zmin* *zmax*]
  temperature          *int*
```

The *density* keyword computes ion and charge density profiles in the z direction with data sampled every *frequency* step. The file is then regularly updated. *num_bins* is the number of bins in the histogram, that starts at *zmin* (by default 0.0) and ends at *zmax* (by default the box length in z).
The *restart* keyword generates data files from which one can restart a simulation. If the keyword *alternate* is specified, restart files are created and alternatively named *restart.1.out* and *restart.2.out* so that there is always one restart saved. If not, a separate restart file *restart.step.out* is created for each step.
The *energies* keyword generates an energy breakdown that by default contains all pertinent informations for the current simulation. However, one can choose the energy terms to be output by specifying one or more keywords from the following list:
```
Total                     # total energy = kinetic + potential energy
Kinetic                   # kinetic energy
Potential                 # potential energy = intramolecular + van_der_Waals + coulomb energy
Intramolecular            # intramolecular energy = bond + angle + dihedral/impropers + steele energy 
Harmonic_bond             # harmonic bond energy
Harmonic_angle            # harmonic angle energy
Dihedrals                 # dihedrals + impropers energy
Coulomb                   # coulomb energy = lr + sr + keq0 + self + molecule energies
Coulomb_lr                # coulomb long range energy
Coulomb_sr                # coulomb short range energy
Coulomb_keq0              # coulomb k=0 energy
Coulomb_self              # coulomb self energy
Coulomb_molecule          # coulomb correction due to molecules
van_der_Waals             # van der Waals energy = LJ + FT + DAIM + XFT energies
Lennard-Jones             # Lennard-Jones atoms contribution to the energy
Lennard-Jones_cutoff      # Lennard-Jones cutoff contribution to the energy
Fumi-Tosi                 # Fumi-Tosi atoms contribution to the energy
Fumi-Tosi_molecule        # Fumi-Tosi correction due to molecules
XFT                       # Extended Fumi-Tosi atoms contribution to the energy
XFT_molecule              # Extended Fumi-Tosi correction due to molecules
DAIM                      # DAIM energy
DAIM_self                 # DAIM self energy
Steele                    # steele energy
Hamiltonian               # hamiltonian (conserved energy) 
Thermostat                # thermostat energy
Barostat                  # barostat energy
Plumed_bias               # bias due to PLUMED
Piston_work               # work made by the pistons
Thomas_Fermi              # Thomas-Fermi additional energy
External_field            # external field additional energy
``` 

The output files are discussed [here](output-files).

   
