lipid_analysis.py

#-------------------------------------------------------------------------
# This file is part of the ost_pymodules project (https://github.com/njohner/ost_pymodules).
#
# Copyright 2015 Niklaus Johner
#
# ost_pymodules is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
#(at your option) any later version.
#
# ost_pymodules is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with ost_pymodules.  If not, see <http://www.gnu.org/licenses/>.
#-------------------------------------------------------------------------
"""
.. codeauthor:: Niklaus Johner <niklaus.johner@a3.epfl.ch>

This module contains functions to determine lipid tilt and splay angles and
Calculate the elastic properties of the membrane from there.

If you use this code please cite ref. [1]_, where the method and its implementation
is described in details or ref. [2]_, where the method is applied to complex lipidic phases.

To ensure proper treatment of the periodic boundary conditions, the trajectory
should first be extended to neighboring unit cells, and then aligned.
The tilts and splays are then calculated only for the central cell, the others
being only used to ascertain correct treatment of the PBCs.

References
-------------
.. [1] Niklaus Johner, D. Harries and G. Khelashvili,
       "Implementation of a methodology for determining elastic properties of lipid assemblies from molecular dynamics simulations",
       submitted to BMC bioinformatics (2015).
.. [2] Niklaus Johner, D. Harries and G. Khelashvili,
       "Curvature and lipid packing modulate the elastic properties of lipid assemblies: comparing the HII and lamellar phases."
       The Journal of Physical Chemistry Letters 5, no. 23 (2014), 4201-6.
"""

try:
    from ost import *
    import time
    import numpy as npy
    import os
    import math
    from scipy.optimize import curve_fit
    import matplotlib as mpl
    try:
        gui.dng
    except:
        mpl.use('Agg')
    import matplotlib.pyplot as plt
    import entity_alg
    import trajectory_utilities
    import surface_alg
    import file_utilities
except ImportError:
    print 'could not import at least one of the modules nedded: ost, time, numpy, os, math, entity_alg,trajectory_utilities,surface_alg,file_utilities'


__all__ = ('GetBoundaryBetweenViews', 'AssignNormalsToLipids', 'AnalyzeLipidTilts',
           'AnalyzeLipidSplays', 'AnalyzeLipidTiltAndSplay',
           'WriteTiltDict', 'WriteSplayDict',
           'FitTiltDistribution', 'FitSplayDistribution', 'AnalyzeAreaPerLipid',
           'ExtractTiltAndSplayModuli')


def _CalculateSplayAngle(v11, v12, v21, v22, v1p, v2p, n1, n2, distance_cutoff):
    """
    Calculates the splay angle for a pair of lipids.

    :param v11: headgroup of the first lipid
    :param v12: terminal tail atoms of the first lipid
    :param v21: headgroup of the second lipid
    :param v22: terminal tail atoms of the second lipid
    :param v1p: Atoms of the first lipid situated at the neutral/pivotal plane
    :param v1p: Atoms of the second lipid situated at the neutral/pivotal plane
    :param distance_cutoff: The maximal distance between lipids considered for splay calculation.

    :type v11: :class:`~ost.mol.EntityView`
    :type v12: :class:`~ost.mol.EntityView`
    :type v21: :class:`~ost.mol.EntityView`
    :type v22: :class:`~ost.mol.EntityView`
    :type v1p: :class:`~ost.mol.EntityView`
    :type v2p: :class:`~ost.mol.EntityView`
    :type distance_cutoff: :class: `float`

    :return: returns a tuple of floats composed of the splay angle and the distance between the lipids
    :rtype: (:class:`float`,:class:`float`)
    """
    if geom.Dot(n1, n2) < 0:
        return None
    x = v2p.GetCenterOfMass() - v1p.GetCenterOfMass()
    d = geom.Length(x)
    x = x - geom.Dot(x, n1) * n1
    if d > distance_cutoff:
        return None
    x = geom.Normalize(x)
    v1 = geom.Normalize(v11.GetCenterOfMass() - v12.GetCenterOfMass())
    v2 = geom.Normalize(v21.GetCenterOfMass() - v22.GetCenterOfMass())
    return ((geom.Dot(v2, x) - geom.Dot(v1, x) - geom.Dot(n2, x) + geom.Dot(n1, x)) / d, d)


CalculateInterfaceFromTraj = trajectory_utilities.CalculateInterfaceFromTraj


def _AssignNormalsFromSurfaceToResidues(t, sele, surface, within_size=10):
    """
    This function assigns a normal vector to each residue for each frame of a trajectory
    from the closest point on the surface.

    :param t: the trajectory
    :param sele: the selection to which normals will be assigned.
    :param surface: the surface. Each atom of the surface should have an associated
                    normal vector as Vec3 property "n".
    :param within_size: Size of surrounding used to find the closest atom on the surface.
                        This parameter only optimizes the speed of the calculation

    :type t: :class:`~ost.mol.CoordGroupHandle`
    :type sele: :class:`~ost.mol.EntityView`
    :type surface: :class:`~ost.mol.EntityView`
    :type within_size: :class:`float`

    :return: A list containing one **Vec3List** for each residue in **sele**.
    """
    surface = surface.Select('')
    ln = [geom.Vec3List() for r in sele.residues]
    res_view_list = [r.Select('') for r in sele.residues]
    for i in range(t.GetFrameCount()):
        t.CopyFrame(i)
        for j, r in enumerate(res_view_list):
            within = surface.FindWithin(r.center_of_atoms, within_size)
            if not len(within) == 0:
                a = entity_alg.FindClosestAtom(
                    r, mol.CreateViewFromAtoms(within))
            else:
                a = entity_alg.FindClosestAtom(r, surface)
            ln[j].append(a.GetVec3Prop('n'))
    return ln


def _CalculateTilts(t, lipids, normals, head_sele, tail_sele,
                    prot_cm=None, bool_prop=''):
    """
    This function calculates the lipid tilts.
    :param t: the trajectory
    :param lipids: the selection for which tilts will be calculated.
    :param head_sele: the selection used for the headgroups.
    :param tail_sele: the selection used for the tails.
    :param prot_cm: A point from which the distance will be calculated for each tilt.
                    This is typically the position of the center of mass of a protein.

    :type t: :class:`~ost.mol.CoordGroupHandle`
    :type lipids: :class:`~ost.mol.EntityView`
    :type head_sele: :class:`str`
    :type tail_sele: :class:`str`
    :type bool_prop: :class:`str`
    :type prot_cm: :class:`~ost.mol.EntityView`

    :return: A tuple of arrays **(tilts, prot_dist)**. Each array has the shape N\ :subscript:`Lipids`\ x N\ :subscript:`Frames`
     If **prot_cm** is **None**, the second list is empty.
    :rtype: (:class:`npy.array`,:class:`npy.array`)
    """
    tilts = []
    prot_dist = []
    for i, r in enumerate(lipids.residues):
        if bool_prop and not r.GetBoolProp(bool_prop):
            continue
        p1 = mol.alg.AnalyzeCenterOfMassPos(t, r.Select(head_sele))
        p2 = mol.alg.AnalyzeCenterOfMassPos(t, r.Select(tail_sele))
        n = normals[i]
        tilts.append([geom.Angle(ni, p1i - p2i)
                      for ni, p1i, p2i in zip(n, p1, p2)])
        if prot_cm:
            prot_dist.append([min(geom.Length(geom.Vec2(cm - el1)), geom.Length(
                geom.Vec2(cm - el2))) for cm, el1, el2 in zip(prot_cm, p1, p2)])
    return (npy.array(tilts), npy.array(prot_dist))


def GetBoundaryBetweenViews(t, waters, lipids, outdir='', density_cutoff=None,
                            stride=1, within_size_normals=5.0, filename_basis=''):
    """
    This function determines the interface between two views, typically the
    lipid-water interface and assigns normals to every point on the surface.

    :param t: the trajectory
    :param waters: First view
    :param lipids: Second view
    :param outdir: Path to output directory
    :param density_cutoff: Interface will not be calculated for regions where
     the density is lower than this cutoff.
    :param stride: stride used to calculate the average density from the trajectory.
    :param within_size_normals: radius of the patch used to determine
     the normals on the water-lipid interface.
    :param filename_basis: used as first part in the name of all the files generated.

    :type t: :class:`~ost.mol.CoordGroupHandle`
    :type waters: :class:`~ost.mol.EntityView`
    :type lipids: :class:`~ost.mol.EntityView`
    :type outdir: :class:`str`
    :type density_cutoff: :class:`float`
    :type stride: :class:`int`
    :type within_size: :class:`float`
    :type filename_basis: :class:`str`

    :return: A tuple **(water_filtered,lipid_filtered,b_eh)** containing the density for
     the first and second views and an entity for the boundary between the two views.
     Every atom in the boundary has an associated normal vector set as a Vec3 property 'n'.

    WARNING: Interface was changed. Taking 2 views instead of a list of lipid
    names and water name. The order of the two views was also inversed.
    I also removed the PBC, cell_center and cell_size parameters

    """

    t0 = time.time()
    eh = t.GetEntity()
    t2 = t.Filter(eh.Select(''), stride=stride)
    print 'Calculating densities on', t2.GetFrameCount(), 'frames'
    (water_filtered, lipid_filtered, boundary_filtered) = CalculateInterfaceFromTraj(
        t2, waters, lipids, False, None, None, 10, density_cutoff, outdir, filename_basis)
    b_eh = entity_alg.CreateEntityFromVec3List(boundary_filtered)
    print 'calculating the normals', time.time() - t0
    surface_alg.CalculateNormals(b_eh, within_size_normals, False, None, None)
    surface_alg.OrientNormalsAlongDensity(b_eh, lipid_filtered)
    if outdir:
        io.SavePDB(b_eh, os.path.join(outdir, filename_basis + 'boundary.pdb'))
        entity_alg.WriteVec3Prop(b_eh, "n", os.path.join(
            outdir, filename_basis + 'boundary_normals.txt'), index=True)
    return (water_filtered, lipid_filtered, b_eh)


def AssignNormalsToLipids(t, eh, b_eh, lipid_names, head_group_dict):
    """
    :param t: The trajectory
    :param eh: The associated entity
    :param b_eh: The surface from which the normals will get assigned to the lipids.
    :param lipid_names: a list of the lipid residue names to which normals will be assigned.
    :param head_group_dict: a dictionary containing the selections defining the headgroup
     for each lipid type. There should be one entry for each residue name in **lipid_names**

    :type t: :class:`~ost.mol.CoordGroupHandle`
    :type eh: :class:`~ost.mol.Entity`
    :type b_eh: :class:`~ost.mol.Entity`
    :type lipid_names: :class:`list`
    :type head_group_dict: :class:`dict`

    :return: A dictionary with one entry for each residue name in **lipid_names**.
     Each element in the dictionary is a :class:`list`\ (:class:`~ost.geom.Vec3List`\ ).
     Each element in the list corresponds to one residue and each
     :class:`~ost.geom.Vec3` in the :class:`~ost.geom.Vec3List` is the normal
     for one frame of the trajectory.
    """
    t0 = time.time()
    # Now we assign a normal for each lipid in each frame
    print 'assigning normals to lipids', lipid_names
    lipid_normal_dict = {}
    for lipid_name in lipid_names:
        print 'rname=' + lipid_name + ' and ' + head_group_dict[lipid_name]
        t0 = time.time()
        lipids = eh.Select('rname=' + lipid_name +
                           ' and ' + head_group_dict[lipid_name])
        print lipids.GetResidueCount()
        lipid_normal_dict[lipid_name] = _AssignNormalsFromSurfaceToResidues(
            t, lipids, b_eh)
        print 'done for', lipid_name, 'in', time.time() - t0, 'seconds'
    return lipid_normal_dict


def AnalyzeLipidTilts(t, eh, lipid_names, lipid_normal_dict, head_group_dict,
                      tail_dict, prot_cm=None, bool_prop=''):
    """
    This function calculates the lipid tilts from a trajectory.

    :param t: The trajectory

    :param eh: The associated entity

    :param lipid_names: List of the residue names of the different lipids in the system

    :param lipid_normal_dict: Dictionary of normal vectors.
     One entry for every lipid type (element in lipid_names)
     Every entry is a :class:`list`\ (:class:`~ost.geom.Vec3List`\ ) of normals
     for every frame for every lipid of that type
     (size of list:N\ :subscript:`Lipids`\ x N\ :subscript:`Frames`).

    :param head_group_dict: Dictionary containing a selection string for each lipid type
     that is used to determine the position of the lipid headgroups (center of mass of the selection).

    :param tail_dict: Dictionary containing a selection string for each lipid type
     that is used to determine the position of the lipid tails (center of mass of the selection).

    :param prot_cm: A list of position (one for each frame).
     If specified the each tilt the distance between this position
     and the lipid in question will also be returned. This is typically used
     to calculate local properties of the membrane around an insertion.

    :param bool_prop: Boolean property assigned to lipids to determine whether
     they should be considered in the tilt calculations.
     This is typically used to treat the periodic boundary conditions,
     to differentiate lipids from the central unit cell, for which tilt and
     splay are calculated, from the lipids from neighboring unit cells,
     used only to ensure correct treatment of PBC.

    :type t: :class:`~ost.mol.CoordGroupHandle`
    :type eh: :class:`~ost.mol.EntityHandle`
    :type lipid_normal_dict: :class:`dict`
    :type lipid_names: :class:`str`
    :type head_group_dict: :class:`dict`
    :type tail_dict: :class:`dict`
    :type prot_cm: :class:`~ost.geom.Vec3List`
    :type bool_prop: :class:`bool`

    :return: Dictionary of lipid tilts. One entry for every lipid type (element in **lipid_names**).
             Every entry is a :class:`list` with two elements. The first one is a :class:`list`\ (:class:`~ost.geom.FloatList`\ ) of tilts
             for every frame for every lipid of that type (size of list:N\ :subscript:`Lipids`\ x N\ :subscript:`Frames`).
             The second element contains the list of distances to **prot_cm** if it was defined and is empty otherwise.

    WARNING: Removed parameters PBC, cell_center, cell_size
    """
    lipid_tilt_dict = {}
    print "analyzing lipid tilts for", lipid_names
    for ln in lipid_names:
        print "starting for", ln
        lipids = eh.Select('rname=' + ln)
        lipid_tilt_dict[ln] = _CalculateTilts(t, lipids, lipid_normal_dict[
                                              ln], head_group_dict[ln], tail_dict[ln], prot_cm, bool_prop)
    return lipid_tilt_dict


def AnalyzeLipidSplays(t, eh, lipid_names, head_group_dict, tail_dict,
                       lipid_normal_dict, lipid_tilt_dict, distance_sele_dict,
                       distance_cutoff=10, bool_prop='', do_prot_dist=False):
    """
    This function calculates the lipid splays from a trajectory.

    :param t: The trajectory
    :param eh: The associated entity
    :param lipid_names: A list of the residue names of the different lipids in the system
    :param head_group_dict: Dictionary containing a selection string for each lipid type
     that is used to determine the position of the lipid headgroups (center of mass of the selection).
    :param tail_dict: Dictionary containing a selection string for each lipid type
     that is used to determine the position of the lipid tails (center of mass of the selection).
    :param lipid_normal_dict: Dictionary of normal vectors. One entry for every lipid type (element in lipid_names)
                              Every entry is a :class:`list`\ (:class:`~ost.geom.Vec3List`\ ) of normals for every frame
                              for every lipid of that type (size of list:N\ :subscript:`Lipids`\ x N\ :subscript:`Frames`).
    :param lipid_tilt_dict: Dictionary of lipid tilts. One entry for every lipid type (element in lipid_names)
                            Every entry is a :class:`list`\ (:class:`~ost.geom.FloatList`\ ) of tilts for every frame
                            for every lipid of that type (size of list:N\ :subscript:`Lipids`\ x N\ :subscript:`Frames`).
    :param distance_sele_dict: Dictionary containing a selection string for each lipid type that is used
     to calculate the distance between lipids (center of mass distance). The center of mass of these selections should lie
     on the neutral plane.
    :param distance_cutoff: Lipid pairs further apart than this distance will not be considered for splay calculation.
    :param bool_prop: Boolean property assigned to lipids to determine whether they should be considered in the splay calculations.
     This is typically used to treat the periodic boundary conditions, to differentiate lipids from the central unit cell, for which tilt and
     splay are calculated, from the lipids from neighboring unit cells, used only to ensure correct treatment of PBC.

    :return: Dictionary of splays, containing one entry for each possible lipid pairs

    WARNING: Changed the name of that function from AnalyzeLipidSplay to AnalyzeLipidSplays.
    """
    t0 = time.time()
    splay_dict = {}

    lipid_sele_string = 'rname=' + \
        ','.join([rname for rname in distance_sele_dict.keys()])
    distance_sele_string = ' or '.join(['(rname={0} and {1})'.format(
        rname, distance_sele_dict[rname]) for rname in distance_sele_dict.keys()])
    print 'selecting all the lipids using', lipid_sele_string
    print 'selection for distance calculation', distance_sele_string
    lipids = eh.Select(lipid_sele_string)
    distance_sele = lipids.Select(distance_sele_string)
    print bool_prop, lipid_names
    head_group_view_list_dict = {}
    tail_view_list_dict = {}
    distance_view_list_dict = {}
    for ln in lipid_names:
        v = lipids.Select('rname=' + ln)
        head_group_view_list_dict[ln] = [
            r.Select(head_group_dict[ln]) for r in v.residues]
        tail_view_list_dict[ln] = [r.Select(tail_dict[ln]) for r in v.residues]
        distance_view_list_dict[ln] = [
            r.Select(distance_sele_dict[ln]) for r in v.residues]
    splay_dict = {}  # We prepare the dictionary that will contain the splays for each type of lipid pair
    for ln1 in lipid_names:
        for ln2 in lipid_names:
            splay_dict[ln1 + '-' + ln2] = []
    for f in range(t.GetFrameCount()):
        if f % 50 == 0:
            print "frame {0} in {1} seconds".format(f, time.time() - t0)
        t.CopyFrame(f)
        for r1 in lipids.residues:
            if not r1.GetBoolProp(bool_prop):
                continue
            ln1 = r1.GetName()
            i = r1.GetIntProp('index')
            within = distance_sele.Select("{0}<>[cname={1} and rnum={2}]".format(
                distance_cutoff, r1.chain.name, r1.number.num))
            within = within.Select("not (cname={0} and rnum={1})".format(
                r1.chain.name, r1.number.num))
            # if within.GetResidueCount()<2:continue
            v11 = head_group_view_list_dict[ln1][i]
            v12 = tail_view_list_dict[ln1][i]
            v13 = distance_view_list_dict[ln1][i]
            n1 = lipid_normal_dict[ln1][i][f]
            a1 = lipid_tilt_dict[ln1][0][i][f]
            for r2 in within.residues:
                ln2 = r2.name
                j = r2.GetIntProp('index')
                v21 = head_group_view_list_dict[ln2][j]
                v22 = tail_view_list_dict[ln2][j]
                v23 = distance_view_list_dict[ln2][j]
                n2 = lipid_normal_dict[ln2][j][f]
                a2 = lipid_tilt_dict[ln2][0][j][f]
                if a1 > 0.8 or a2 > 0.8:
                    continue
                s = _CalculateSplayAngle(
                    v11, v12, v21, v22, v13, v23, n1, n2, distance_cutoff)
                # Distance to protein
                if do_prot_dist:
                    s += (0.5*(lipid_tilt_dict[ln1][0][i][f]+lipid_tilt_dict[ln2][1][j][f]),)
                if s:
                    splay_dict[ln1 + '-' + ln2].append(s)
    # We finalize the dictionary by merging equivalent entries
    nsplays = 0
    splay_dict_def = {}
    for ln1 in lipid_names:
        for ln2 in lipid_names:
            key = ln1 + '-' + ln2
            key2 = ln2 + '-' + ln1
            nsplays += len(splay_dict[key])
            if key2 in splay_dict_def.keys():
                splay_dict_def[key2].extend(splay_dict[key])
            else:
                splay_dict_def[key] = splay_dict[key]
    print 'done in {0} seconds. Computed {1} splays'.format(time.time() - t0, nsplays)
    return splay_dict_def


def AnalyzeLipidTiltAndSplay(t, lipid_names, head_group_dict, tail_dict,
                             distance_cutoff=10.0, within_size_normals=10.0,
                             distance_sele_dict=None, water_name='TIP3', outdir='',
                             density_cutoff=None, prot_sele=None, density_stride=10,
                             tilt_bool_prop='', splay_bool_prop='', filename_basis='',
                             sele_dict=None):
    """
    This function is a wrapper to determine the membrane elastic moduli from the lipid titls and splays.
    Periodic boundary conditions are not treated explicitely here and should be treated as suggested in the
    description of this module.

    :param t: The trajectory
    :param lipid_names: List of the residue names of the different lipids in the system
    :param head_group_dict: Dictionary containing a selection string for each lipid type
     that is used to determine the position of the lipid headgroups (center of mass of the selection).
    :param tail_dict: Dictionary containing a selection string for each lipid type
     that is used to determine the position of the lipid tails (center of mass of the selection).
    :param distance_cutoff: Lipid pairs further apart than this distance will not be considered for splay calculation.
    :param within_size_normals: radius of the patch used to determine the normals on the water-lipid interface.
    :param distance_sele_dict: Dictionary containing a selection string for each lipid type that is used
     to calculate the distance between lipids (center of mass distance). The center of mass of these selections should lie
     on the neutral plane.
    :param water_name: Residue name of the waters (used to calculate the water-lipid interface).
    :param outdir: Path to output directory. If none, no files will be written.
    :param density_cutoff: Interface will not be calculated for regions where the density is lower than this cutoff.
    :param prot_sele: Selection string used to determine the position of the protein. This is used to calculate the distance
     between lipids and the protein which will be returned together with splays and tilts and can be used to determine the vatiation
     in membrane properties around a protein
    :param density_stride: Stride to be used for the calculation of the average lipid and water densities used to determine the interface.
     Using every single frame can slow down the calculation.
    :param tilt_bool_prop: Boolean property assigned to lipids to determine whether they should be considered in the titl calculations.
     This is typically used to treat the periodic boundary conditions, to differentiate lipids from the central unit cell, for which tilt and
     splay are calculated, from the lipids from neighboring unit cells, used only to ensure correct treatment of PBC.
    :param splay_bool_prop: Same as **tilt_bool_prop** but for the calculation of splays
    :param filename_basis: used as first part in the name of all the files generated.
    :param sele_dict: Dictionary containing selection strings used to separate the system in several parts for the calculation
     This can be used for example to make the tilt and splay calculations separately for each leaflet of a bilayer.
     In such a case **sele_dict** would be something like **sele_dict={"upper":"z>0","lower":"z<0"}**

    :type t: :class:`~ost.mol.CoordGroupHandle`
    :type lipid_names: :class:`str`
    :type head_group_dict: :class:`dict`
    :type tail_dict: :class:`dict`
    :type distance_cutoff: :class:`float`
    :type within_size_normals: :class:`float`
    :type distance_sele_dict: :class:`dict`
    :type outdir: :class:`str`
    :type density_cutoff: :class:`float`
    :type prot_sele: :class:`str`
    :type density_stride: :class:`int`
    :type tilt_bool_prop: :class:`bool`
    :type splay_bool_prop: :class:`bool`
    :type filename_basis: :class:`str`
    :type sele_dict: :class:`dict`

    :return: A tuple **(lipid_tilt_dict,lipid_normal_dict,splay_dict,b_eh)**, where **lipid_tilt_dict**,
            **lipid_normal_dict** and **lipid_splay_dict** are dictionaries with keys corresponding to the elements
            in **sele_dict**. For more information about **lipid_tilt_dict** and **lipid_splay_dict**, refer to the
            the documentation for the **AnalyzeLipidTilts** and **AnalyzeLipidSplays** functions.

    WARNING: Removed parameters PBC, cell_center, cell_size
    """
    if distance_sele_dict is None:
        distance_sele_dict = {}
    if sele_dict is None:
        sele_dict = {}

    t0 = time.time()
    eh = t.GetEntity()
    t.CopyFrame(0)
    # first we build the interafce between water and membrane
    # and assign a normal for each lipid in each frame
    print 'Generating the water and lipid densities and boundary surface'
    lipid_sele = 'rname=' + ','.join(lipid_names)
    water_sele = 'rname=' + water_name
    lipids = eh.Select(lipid_sele)
    waters = eh.Select(water_sele)
    (water_filtered, lipid_filtered, b_eh) = GetBoundaryBetweenViews(t, waters, lipids, outdir, density_cutoff,
                                                                     density_stride, within_size_normals, filename_basis)
    # Protein center of mass
    if prot_sele:
        prot_cm = mol.alg.AnalyzeCenterOfMassPos(t, eh.Select(prot_sele))
        do_prot_dist = True
    else:
        prot_cm = None
        do_prot_dist = False
    # From here we separate into the sub categories defined in sele_dict,
    # which could be upper and lower leaflets
    s = ['(rname=' + ln + ' and ' + head_group_dict[ln] + ')' for ln in lipid_names]
    lipid_sele_dict = {}
    if not sele_dict:
        sele_dict['all'] = ''
    for sele_name in sele_dict:
        if sele_dict[sele_name]:
            lipid_sele_dict[sele_name] = ' or '.join(
                s) + ' and ' + sele_dict[sele_name]
        else:
            lipid_sele_dict[sele_name] = ' or '.join(s)
    print sele_dict, lipid_sele_dict
    lipid_normal_dict = {}
    lipid_tilt_dict = {}
    if tilt_bool_prop:
        lipid_tilt_dict_sele = {}
    splay_dict = {}
    for sele_name in sele_dict:
        lipid_sele = lipid_sele_dict[sele_name]
        sele = sele_dict[sele_name]
        sele_ev = eh.Select(lipid_sele, mol.MATCH_RESIDUES)
        for ln in lipid_names:
            v = sele_ev.Select('rname={0}'.format(ln))
            for i, r in enumerate(v.residues):
                r.SetIntProp('index', i)
        print 'Assigning normals for', sele_name, sele, lipid_sele
        lipid_normal_dict[sele_name] = AssignNormalsToLipids(
            t, sele_ev, b_eh.Select(sele), lipid_names, head_group_dict)
        print 'lipid normal dict', lipid_normal_dict[sele_name].keys()
        print 'Done in', time.time() - t0, 'seconds'
        t0 = time.time()
        print 'calculating tilts for', sele
        lipid_tilt_dict[sele_name] = AnalyzeLipidTilts(t, eh.Select(lipid_sele, mol.MATCH_RESIDUES), lipid_names, lipid_normal_dict[sele_name],
                                                       head_group_dict, tail_dict, prot_cm)
        if tilt_bool_prop:
            print "Caalculating tilts again only for the lipids with bool prop", tilt_bool_prop
            lipid_tilt_dict_sele[sele_name] = AnalyzeLipidTilts(t, eh.Select(lipid_sele, mol.MATCH_RESIDUES), lipid_names, lipid_normal_dict[sele_name],
                                                                head_group_dict, tail_dict, prot_cm, tilt_bool_prop)
        print 'Done in', time.time() - t0, 'seconds'
        t0 = time.time()
        print 'calculating splay for', sele
        splay_dict[sele_name] = AnalyzeLipidSplays(t, eh.Select(lipid_sele, mol.MATCH_RESIDUES), lipid_names, head_group_dict, tail_dict,
                                                   lipid_normal_dict[sele_name], lipid_tilt_dict[sele_name], distance_sele_dict, distance_cutoff, splay_bool_prop, do_prot_dist)
        print 'Done in', time.time() - t0, 'seconds'
        if outdir:
            try:
                WriteSplayDict(
                    splay_dict[sele_name], outdir, filename_basis + sele_name + '_', do_prot_dist)
            except IOError:
                print 'could not write dict'
            if tilt_bool_prop:
                WriteTiltDict(lipid_tilt_dict_sele[
                              sele_name], outdir, filename_basis + sele_name + '_')
            else:
                WriteTiltDict(lipid_tilt_dict[
                              sele_name], outdir, filename_basis + sele_name + '_')
    if tilt_bool_prop:
        return (lipid_tilt_dict_sele, lipid_normal_dict, splay_dict, b_eh)
    else:
        return (lipid_tilt_dict, lipid_normal_dict, splay_dict, b_eh)


def WriteTiltDict(lipid_tilt_dict, outdir, filename_basis=''):
    """
    This function writes out the *lipid_tilt_dict* to the firectory *outdir*.
    It writes out one file for every key in *lipid_tilt_dict*, i.e. for every lipid type.
    Filenames are preceded by the *filename_basis*.

    :param lipid_tilt_dict: the lipid tilt dictionary as produced by :func:`AnalyzeLipidTiltAndSplay`.
    :param outdir: the directory to which the files will be written
    :param filename_basis: Will be prepended to all file names.

    :type lipid_tilt_dict: :class:`dict`
    :type outdir: :class:`str`
    :type filename_basis: :class:`str`
    """
    for ln in lipid_tilt_dict:
        tl = [180. * el / math.pi for lt in lipid_tilt_dict[ln][0]
              for el in lt]
        dl = [el for lt in lipid_tilt_dict[ln][1] for el in lt]
        if len(dl) != 0:
            if len(dl) != len(tl):
                print 'not the same number of tilts and distances', len(tl), len(dl)
            ll = [tl, dl]
            titles = ['tilt', 'dist from prot']
        else:
            ll = [tl]
            titles = ['tilt']
        file_utilities.WriteListOfListsInColumns(titles, ll, os.path.join(
            outdir, filename_basis + 'tilt_' + ln + '.txt'))


def WriteSplayDict(splay_dict, outdir, filename_basis, prot_dist_flag):
    """
    This function writes out the *splay_dict* to the firectory *outdir*.
    It writes out one file for every key in *splay_dict*, i.e. for every lipid type pair.
    Filenames are preceded by the *filename_basis*.

    :param splay_dict: the lipid splay dictionary as produced by :func:`AnalyzeLipidTiltAndSplay`.
    :param outdir: the directory to which the files will be written
    :param filename_basis: Will be prepended to all file names.

    :type splay_dict: :class:`dict`
    :type outdir: :class:`str`
    :type filename_basis: :class:`str`
    """
    if len(splay_dict.keys()) == 0:
        return
    for ln in splay_dict:
        if len(splay_dict[ln]) == 0:
            continue
        if prot_dist_flag:
            sl = [[el[0], el[1], el[2]] for el in splay_dict[ln]]
            file_utilities.WriteListOfListsInLines(['splay', 'dist', 'ProtDist'], sl, os.path.join(
                outdir, filename_basis + 'splay_' + ln + '.txt'))
        else:
            sl = [[el[0], el[1]] for el in splay_dict[ln]]
            file_utilities.WriteListOfListsInLines(['splay', 'dist'], sl, os.path.join(
                outdir, filename_basis + 'splay_' + ln + '.txt'))
    return


def _gauss(x, *p):
    A, mu, sigma = p
    return A * npy.exp(-(x - mu)**2 / (2. * sigma**2))


def _parabole(x, *p):
    a, b, x0 = p
    return a + b * (x - x0)**2.0


def _centered_parabole(x, *p):
    a, b = p
    return a + b * x**2.0


def _FindIndexOfClosestValue(l, v):
    return min(enumerate(l), key=lambda x: abs(x[1] - v))[0]


def _FitGaussian(bincenters, pa):
    mu0 = npy.sum(bincenters * pa) / npy.sum(pa)
    A0 = max(pa)
    sigma0 = npy.sqrt(npy.sum(((bincenters - mu0)**2.0) * pa) / npy.sum(pa))
    (A, mu, sigma), v = curve_fit(_gauss, bincenters, pa, [A0, mu0, sigma0])
    return A, mu, abs(sigma)


def _PlotGaussian(bincenters, pa, A, mu, sigma, outfile, title='', xlabel=''):
    plt.figure()
    plt.plot(bincenters, pa, 'o')
    plt.plot(bincenters, _gauss(bincenters, A, mu, sigma), '-', color='g',
             label='$\mu={0}\ ;\ \sigma={1}$'.format(round(mu, int(1 / mu) + 1), round(sigma, int(1 / sigma) + 1)))
    plt.vlines(mu, 0, plt.ylim()[1], linestyle='--', color='r')
    plt.vlines([mu + sigma, mu - sigma], 0, plt.ylim()
               [1], linestyle='--', color='c')
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel('Probability')
    plt.savefig(outfile)
    plt.close()


def _FitParabole(bincenters, fa, fitting_range, centered=False):
    first = _FindIndexOfClosestValue(bincenters, fitting_range[0])
    last = _FindIndexOfClosestValue(bincenters, fitting_range[1])
    mask = fa != npy.inf
    a = min(fa)
    x0 = bincenters[npy.argmin(fa)]
    xm = bincenters[mask][npy.argmax(fa[mask])]
    fm = max(fa[mask])
    b = (fm - a) / (xm - x0)**2.0
    if centered:
        r, v = curve_fit(_centered_parabole, bincenters[
                         first:last], fa[first:last], [a, b])
    else:
        r, v = curve_fit(_parabole, bincenters[first:last], fa[
                         first:last], [a, b, x0])
    return r


def _PlotParabola(bincenters, fa, a, b, x0, fitting_range, outfile, title='', xlabel='', ylabel='', fit_label=''):
    plt.figure()
    plt.plot(bincenters, fa, 'o')
    ymax = plt.ylim()[1]
    ymin = plt.ylim()[0]
    plt.plot(bincenters, [_parabole(xi, a, b, x0)
                          for xi in bincenters], '--', color='r', label=fit_label)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.title(title)
    if fit_label:
        plt.legend(loc='best')
    plt.vlines(fitting_range, ymin, ymax, linestyle='--', color='c')
    plt.ylim(ymin, ymax)
    plt.savefig(outfile)
    plt.close()
    return


def ExtractTiltAndSplayModuli(tilt_dict, splay_dict, lipid_area, outdir, nbins=100):
    """
    This function extracts the tilt and splay moduli from the dictionaries of tilts and splays
    obtained from the *AnalyzeLipidTiltAndSplay* function.
    It will first fit a gaussian y=A exp[(x-mu)/sigma^2] to the distribution of tilts and the distribution of splays
    to determine the fitting range then used to fit the corresponding potential of mean force (PMF).
    Different fitting ranges are used to estimate the error on the extracted tilt and splay moduli.
    The function will calculate one tilt modulus for each lipid species and one splay modulus for each pair
    of lipid species. It will then combine these to calculate the overall tilt modulus and splay modulus (bending rigidity).
    More details about this procedure can be found in ref. [2]_
    The function returns a dictionary containing the tilt and splay moduli with their uncertainties as tuples (k,delta_k).

    :param tilt_dict: A dictionary of lipid tilts as returned by the *AnalyzeLipidTiltAndSplay* function.
    :param splay_dict: A dictionary of lipid splays as returned by the *AnalyzeLipidTiltAndSplay* function.
    :param lipid_area: The area per lipid.
    :param outdir:  the directory to which output files will be written, i.e. plots of tilt and splay distributions and PMFs.
                    and text files containing the elastic moduli
    :param nbins:    The number of bins used when determining the distributions of tilts and splays

    :type tilt_dict: :class:`dict`
    :type splay_dict: :class:`dict`
    :type outdir: :class:`str`
    :type nbins: :class:`int`

    :return: a dictionary containing the tilt and splay moduli with their uncertainties as tuples (k,delta_k).
    :rtype: :class:`dict`

    """
    kdict = {"tilt": {}, "splay": {}}
    for sele_name in tilt_dict.keys():
        kdict["tilt"][sele_name] = {}
        outfile = open(os.path.join(
            outdir, sele_name + '_tilt_constants.txt'), 'w')
        k_list = FloatList()
        deltak_list = FloatList()
        nl_list = FloatList()
        for lipid_name in tilt_dict[sele_name].keys():
            fname = "_".join([sele_name, lipid_name])
            tilt_list = FloatList()
            for el in tilt_dict[sele_name][lipid_name][0]:
                tilt_list.extend(el)
            k, dk, kl = FitTiltDistribution(
                tilt_list, nbins, outdir=outdir, filename_basis=fname, title_complement='for ' + sele_name + " " + lipid_name)
            print "Tilt modulus for {0} {1} is k={2:1.2f} +/- {3:1.2f}.".format(sele_name, lipid_name, k, dk)
            outfile.write(
                ' '.join([lipid_name, str(round(kl[0], 1)), str(round(dk, 1))]) + '\n')
            k_list.append(k)
            deltak_list.append(dk)
            nl_list.append(len(tilt_list))
            kdict["tilt"][sele_name][lipid_name] = (k, dk)
        ntot = npy.sum(nl_list)
        k = npy.sum([ni / ki for ni, ki in zip(nl_list, k_list)]) / ntot
        k = 1. / k
        dk = npy.sum([((k**2.0) / (ki**2.0) * ni / ntot)**2.0 * (dki**2.0)
                      for ni, ki, dki in zip(nl_list, k_list, deltak_list)])
        dk = npy.sqrt(dk)
        outfile.write(
            ' '.join(['agregated', str(round(k, 1)), str(round(dk, 1))]) + '\n')
        outfile.close()
        kdict["tilt"][sele_name]['agregated'] = (k, dk)
        for i in range(2 * len(k_list)):
            plt.close()

    for sele_name in splay_dict.keys():
        kdict["splay"][sele_name] = {}
        outfile = open(os.path.join(
            outdir, sele_name + '_splay_constants.txt'), 'w')
        k_list = FloatList()
        nl_list = FloatList()
        deltak_list = FloatList()
        for key in splay_dict[sele_name]:
            fname = "_".join([sele_name, key])
            splay_list = [el[0] for el in splay_dict[sele_name][key]]
            k, dk, kl = FitSplayDistribution(splay_list, lipid_area, nbins, outdir=outdir,
                                             filename_basis=fname, title_complement='for ' + sele_name + " " + key)
            print "Splay modulus for {0} {1} is k={2:1.2f} +/- {3:1.2f}.".format(sele_name, key, k, dk)
            outfile.write(
                ' '.join([key, str(round(kl[0], 1)), str(round(dk, 1))]) + '\n')
            k_list.append(k)
            deltak_list.append(dk)
            nl_list.append(len(splay_list))
            kdict["splay"][sele_name][key] = (k, dk)
        ntot = npy.sum(nl_list)
        k = npy.sum([ni / ki for ni, ki in zip(nl_list, k_list)]) / ntot
        k = 1. / k
        dk = npy.sum([((k**2.0) / (ki**2.0) * ni / ntot)**2.0 * (dki**2.0)
                      for ni, ki, dki in zip(nl_list, k_list, deltak_list)])
        dk = npy.sqrt(dk)
        outfile.write(
            ' '.join(['agregated', str(round(k, 1)), str(round(dk, 1))]) + '\n')
        outfile.close()
        kdict["splay"][sele_name]['agregated'] = (k, dk)
        for i in range(2 * len(k_list)):
            plt.close()
        return kdict


def FitSplayDistribution(splay_list, lipid_area, nbins=100, x_range=None, outdir='', filename_basis='', title_complement=''):
    """
    This function extracts the bending modulus from a list of splays.
    It will first fit a gaussian y=A exp[(x-mu)/sigma^2] to the distribution of splays to determine
    the fitting range used to fit the potential of mean force (PMF).
    Different fitting ranges are used to estimate the error on the extracted bending rigidity

    :param splay_list: A list of lipid splays.
    :param lipid_area: The area per lipid.
    :param nbins:   The number of bins used when determining the distribution of splays
    :param x_range:  The range in which the distribution will be calculated. Defaults to [mean-3*std,mean+3*std].
    :param outdir:  the directory to which output files will be written, i.e. plots of splay distribution and PMF.
                    if it is not defined, plots will not be generated
    :param filename_basis: Will be prepended to all file names.
    :param title_complement: will be added in the title of the plots

    :type splay_list: :class:`list`
    :type lipid_area: :class:`float`
    :type nbins: :class:`int`
    :type x_range: :class:`tuple` of 2 floats
    :type outdir: :class:`str`
    :type filename_basis: :class:`str`
    :type title_complement: :class:`str`

    :return: A tuple **(K, DeltaK,K_list)**, containing the bending rigidity *K*,
             the estimated uncertainty on *K*, and the list of *K* values obtained from the different fitting ranges.
    :rtype: (:class:`float`,:class:`float`,:class:`list`)
    """
    if not x_range:
        w = npy.std(splay_list)
        m = npy.average(splay_list)
        x_range = [m - 3 * w, m + 3 * w]
    pa = npy.histogram(splay_list, nbins, range=x_range, density=True)
    bincenters = 0.5 * (pa[1][1:] + pa[1][:-1])
    fa = -npy.log(pa[0])
    A, mu, sigma = _FitGaussian(bincenters, pa[0])
    ranges = [(mu - i * sigma, mu + i * sigma)
              for i in [1, 1.25, 1.5, 1.75, 2]]
    print "Gaussian y=A exp[(x-mu)/sigma^2] fitted to tilt distribution with A={0},mu={1} and sigma={2}".format(round(A, 2), round(mu, 2), round(sigma, 2))
    print "Using the following ranges to fit the PMF:", ranges
    res_list = []
    for fitting_range in ranges:
        try:
            r = _FitParabole(bincenters, fa, fitting_range)
        except:
            r = [0, 0]
        res_list.append(r)
    K_list = [2. * res[1] / lipid_area for res in res_list]
    DeltaK = npy.std([el for el in K_list if not el == 0.0])
    K = K_list[0]
    if outdir:
        file_utilities.WriteListOfListsInColumns(['bin', 'distribution', 'pmf'], [bincenters, pa[0], fa], os.path.join(
            outdir, '_'.join([filename_basis, 'splay', 'distribution']) + '.txt'), separator=' ')
        title = 'Splay distribution {0}'.format(title_complement)
        outfile = os.path.join(outdir, '_'.join(
            [filename_basis, 'splay', 'distribution']) + '.png')
        _PlotGaussian(bincenters, pa[0], A, mu, sigma,
                      outfile, title=title, xlabel='Splay')
        outfile = os.path.join(outdir, '_'.join(
            [filename_basis, 'splay', 'fit']) + '.png')
        a, b, x0 = res_list[0]
        title = 'Splay PMF {0}'.format(title_complement)
        _PlotParabola(bincenters, fa, a, b, x0, ranges[0], outfile, title, 'Splay', r'$-\log\left[P(\alpha)\right]$',
                      '$\chi^{{12}}={0}\pm {1}  k_BT$'.format(round(K, int(-math.log10(K)) + 1), round(DeltaK, int(-math.log10(DeltaK)) + 1)))
    return K, DeltaK, K_list


def FitTiltDistribution(tilt_list, nbins=90, x_range=None, outdir='', filename_basis='', title_complement='', degrees=False):
    """
    This function extracts the tilt modulus from a list of lipid tilts.
    It will first fit a gaussian y=A exp[(x-mu)/sigma^2] to the distribution of tilts to determine
    the fitting range used to fit the potential of mean force (PMF).
    Different fitting ranges are used to estimate the error on the extracted tilt modulus

    :param tilt_list: A list of lipid splays.
    :param nbins:   The number of bins used when determining the distribution of splays
    :param x_range:  The range in which the distribution will be calculated. Defaults to [mean-3*std,mean+3*std].
    :param outdir:  the directory to which output files will be written, i.e. plots of splay distribution and PMF.
                    if it is not defined, plots will not be generated
    :param filename_basis: Will be prepended to all file names.
    :param title_complement: will be added in the title of the plots
    :param degrees: Whether plots should be in degrees or radians.

    :type tilt_list: :class:`list`
    :type nbins: :class:`int`
    :type x_range: :class:`tuple` of 2 floats
    :type outdir: :class:`str`
    :type filename_basis: :class:`str`
    :type title_complement: :class:`str`
    :type degrees: :class:`bool`

    :return: A tuple **(Chi, DeltaChi,Chi_list)**, containing the tilt modulus *Chi*,
             the estimated uncertainty on *Chi*, and the list of *Chi* values obtained from the different fitting ranges.
    :rtype: (:class:`float`,:class:`float`,:class:`list`)
    """
    if degrees:
        ac = math.pi / 180.
    else:
        ac = 1
    if x_range is None:
        x_range = [0, math.pi / (2. * ac)]
    pa = npy.histogram(tilt_list, nbins, range=x_range, density=True)
    bincenters = 0.5 * (pa[1][1:] + pa[1][:-1])
    pa = pa[0]
    pa2 = pa / npy.sin(bincenters * ac)
    fa = -npy.log(pa2)
    A, mu, sigma = _FitGaussian(bincenters, pa)
    ranges = [(max(mu - i * sigma, 0), mu + i * sigma)
              for i in [1, 1.25, 1.5, 1.75, 2.0]]
    print "Gaussian y=A exp[(x-mu)/sigma^2] fitted to tilt distribution with A={0},mu={1} and sigma={2}".format(round(A, 2), round(mu, 2), round(sigma, 2))
    print "Using the following ranges to fit the PMF:", ranges
    res_list = [_FitParabole(bincenters, fa, fitting_range, centered=True)
                for fitting_range in ranges]
    K_list = [2. * r[1] / (ac * ac) for r in res_list]
    DeltaK = npy.std(K_list)
    K = K_list[0]
    if outdir:
        file_utilities.WriteListOfListsInColumns(['bin', 'distribution', 'pmf'], [bincenters, pa, fa], os.path.join(
            outdir, '_'.join([filename_basis, 'tilt', 'distribution']) + '.txt'), separator=' ')
        title = 'Tilt distribution {0}'.format(title_complement)
        outfile = os.path.join(outdir, '_'.join(
            [filename_basis, 'tilt', 'distribution']) + '.png')
        _PlotGaussian(bincenters, pa, A, mu, sigma,
                      outfile, title=title, xlabel='Tilt')
        outfile = os.path.join(outdir, '_'.join(
            [filename_basis, 'tilt', 'fit']) + '.png')
        a, b = res_list[0]
        r = [el for el in ranges[0]]
        title = 'Tilt PMF {0}'.format(title_complement)
        _PlotParabola(bincenters, fa, a, b, 0.0, r, outfile, title, 'Tilt', r'$-\log\left[\frac{P(\alpha)}{\sin(\alpha)}\right]$', '$\chi={0}\pm {1} k_BT$'.format(
            round(K, int(-math.log10(K)) + 1), round(DeltaK, int(-math.log10(DeltaK)) + 1)))
    return K, DeltaK, K_list


def AnalyzeAreaPerLipid(t, lipids):
    """
    This function calculates the area per lipid simply from the number of lipids
    and the size of the simulation box.
    The area per lipid is calculated for each frame in the simulation and the average
    is returned. **This is only suitable for bilayers**

    :param t: The trajectory
    :param lipids: Selection of the lipids in the system

    :type t: :class:`~ost.mol.CoordGroupHandle`
    :type lipid: :class:`~ost.mol.EntityView`

    :return: The area per lipid
    :rtype: :class:`float`
    """
    n = lipids.GetResidueCount()
    Al = FloatList()
    for i in range(t.GetFrameCount()):
        f = t.GetFrame(i)
        cell_size = f.GetCellSize()
        Al.append(cell_size[0] * cell_size[1])
    return 2. * npy.average(Al) / float(n)