Asymmetric PC-PE bilayers
=================================================================

This page details the published example (available under /Examples in the git
repositiory). It uses the hoobas package to build the topology, whose documentation
can be found at hoobas.readthedocs.io, and we will concentrate on molecule-wise
semi-grand canonical features here.

Force-field descriptions
-------------------------

Since we perform exchange of angle potentials, we need to create MC type forces for
angles and disable them, on top of the normal angle potentials:

.. code-block:: python

    angle = hoomd.md.angle.cosinesq()
    sgc_angle_saturated = plugin_muMC.angle.cosinesq(name='sat')
    sgc_angle_unsaturated = plugin_muMC.angle.cosinesq(name='unsat')
    for angle_type in builder.angle_types:
        angle.angle_coeff.set(angle_type.typename, k=angle_type['energy_constant'], t0=angle_type['angle_constant'])
        sgc_angle_saturated.angle_coeff.set(angle_type.typename, k=25.0, t0=math.radians(180.0))
        sgc_angle_unsaturated.angle_coeff.set(angle_type.typename, k=45.0, t0=math.radians(120.0))
    sgc_angle_saturated.disable(log=True)
    sgc_angle_unsaturated.disable(log=True)

These forces will be evaluated on MC steps only to compute changes in energies. Similarly,
for every beadtype that can be swapped to (C1, C3, ghost for saturation, Qd, Q0 for headgroup),
we need to define a MC force and disable it.

.. code-block:: python

    sgc_lj_unsaturated = plugin_muMC.pair.lj_mc(nlist=nlist, r_cut=2.5 * 0.62, override='C3', name='sat')
    sgc_lj_unsaturated.disable(log=True)


Molecule ensembles and SGC updaters
------------------------------------

In the input file, we have created a group that contains at least one bead per phospholipid molecule. The
molecule parser will traverse all bonds (and constraints) to construct a full molecule structure. The
resulting structure is a molecule object. This object is then used to hold hash values. To swap chemical
bead nature, we need to describe available states, using a StateDescriptor.

For lipid saturation changes in the MARTINI force-field, this is bead types C1 and C3 and corresponding angle names.
Particle type names and associated MC forces need to be passed in.
The descriptor is then added to build a stateSwitcher object that will perform MC updates.
Note that if all updated types of this simulations (5) are passed in, the hash estimates ceil(log2(NTypes)) * NBeads bits for the hash size.
The maximum value of the hash would therefore be 3 * 8 = 24 bits.
The associated memory footprint to hold chemical potential values would then be 2^24 * 4 bytes.
Memory usage can therefore be minimized by creating multiple updaters if the types are disjoined.
For instance lipid head group is either Q0 or Qd, but will never be C1 or C3 (inversely for acyl tails)


.. code-block:: python

    molecules = plugin_muMC.collectives.Molecules(phospholipids)
    sn1_states = plugin_muMC.update.StateDescriptor(types=['C3', 'C1'], type_force=[sgc_lj_unsaturated, sgc_lj_saturated], angles=[atype120, atype180], angle_force=[sgc_angle_unsaturated, sgc_angle_saturated])
    switcher_sn1 = plugin_muMC.update.ManyStatesSwitcher(group_sn1, molecules, sn1_states, phase=0, name='sn1', stateNames=['unsat', 'sat'])
    switcher_sn1.setRCut(2.5 * 0.62)


Setting chemical potential values
-----------------------------------

The molecule class provides an interface to supplied hash construction in order to translate the binary hash values
to graspable chemistry.
In the input file, we have separated acyl tails updaters, and we therefore obtain different fragments for each tail
The interface provides an iterable, with dict-like notation to obtain the chemical information.
We have named the sn1 updater 'sn1', and the associated name is therefore 'SGCsn1'.
When iterating through the fragments, we can grab the sn1 configuration by element['SGCsn1'].
We named the unsaturations 'unsat' and we can therefore directly count the number of unsaturations by summing all 'unsat's.
Associated hash is given by 'hash' key.


.. code-block:: python

    for element in molecules.fragments():
        fragmentSN1 = element['SGCsn1']
        number_of_unsaturations_in_sn1 = sum([1 for x in fragmentSN1 if x=='unsat'])
        myHash = element['hash']