added custom model files

examples/DSPC_custom_XY/custom_XY_DSPC.py

@@ -0,0 +1,137 @@
+import math
+
+import numpy as np
+
+
+def custom_XY_DSPC(params, bulk_in, bulk_out, contrast):
+    """This function makes a model of a supported DSPC bilayer using volume restricted distribution functions."""
+    # Split up the parameters
+    subRough = params[0]
+    oxideThick = params[1]
+    oxideHydration = params[2]
+    waterThick = params[3]
+    lipidAPM = params[4]
+    lipidCoverage = params[5]
+    bilayerRough = params[6]
+
+    # We are going to need our Neutron scattering cross-sections, plus the component volumes
+    # (the volumes are taken from Armen et al as usual).
+    # Define these first
+
+    # define all the neutron b's.
+    bc = 0.6646e-4  # Carbon
+    bo = 0.5843e-4  # Oxygen
+    bh = -0.3739e-4  # Hydrogen
+    bp = 0.513e-4  # Phosphorus
+    bn = 0.936e-4  # Nitrogen
+
+    # Now make the lipid groups
+    COO = (4 * bo) + (2 * bc)
+    GLYC = (3 * bc) + (5 * bh)
+    CH3 = (2 * bc) + (6 * bh)
+    PO4 = (1 * bp) + (4 * bo)
+    CH2 = (1 * bc) + (2 * bh)
+    CHOL = (5 * bc) + (12 * bh) + (1 * bn)
+
+    # Group these into heads and tails
+    heads = CHOL + PO4 + GLYC + COO
+    tails = (34 * CH2) + (2 * CH3)
+
+    # We need volumes for each. Use literature values
+    vHead = 319
+    vTail = 782
+
+    # Start making our sections. For each we are using a roughened Heaviside to describe our groups.
+    # We will make these as Volume Fractions (i.e. with a height of 1 for full occupancy),
+    # which we will correct later for hydration.
+
+    # Make an array of z values for our model
+    z = np.arange(0, 141)
+
+    # Make our Silicon substrate
+    vfSilicon, siSurf = layer(z, -25, 50, 1, subRough, subRough)
+
+    # Add the Oxide
+    vfOxide, oxSurface = layer(z, siSurf, oxideThick, 1, subRough, subRough)
+
+    # We fill in the water at the end, but there may be a hydration layer between the bilayer and the oxide,
+    # so we start the bilayer stack an appropriate distance away
+    watSurface = oxSurface + waterThick
+
+    # Now make the first lipid head group
+    # Work out the thickness
+    headThick = vHead / lipidAPM
+
+    # ... and make a box for the volume fraction (1 for now, we correct for coverage later)
+    vfHeadL, headLSurface = layer(z, watSurface, headThick, 1, bilayerRough, bilayerRough)
+
+    # ... also do the same for the tails
+    # We'll make both together, so the thickness will be twice the volume
+    tailsThick = (2 * vTail) / lipidAPM
+    vfTails, tailsSurf = layer(z, headLSurface, tailsThick, 1, bilayerRough, bilayerRough)
+
+    # Finally the upper head ...
+    vfHeadR, headSurface = layer(z, tailsSurf, headThick, 1, bilayerRough, bilayerRough)
+
+    # Making the model
+    # We've created the volume fraction profiles corresponding to each of the groups.
+    # We now convert them to SLDs, taking in account of the hydrations to scale the volume fractions
+
+    # 1. Oxide ...
+    vfOxide = vfOxide * (1 - oxideHydration)
+
+    # 2. Lipid ...
+    # Scale both the heads and tails according to overall coverage
+    vfTails = vfTails * lipidCoverage
+    vfHeadL = vfHeadL * lipidCoverage
+    vfHeadR = vfHeadR * lipidCoverage
+
+    # Some extra work to deal with head hydration, which we take to be an additional 30% always
+    vfHeadL = vfHeadL * 0.7
+    vfHeadR = vfHeadR * 0.7
+
+    # Make a total Volume Fraction across the whole interface
+    vfTot = vfSilicon + vfOxide + vfHeadL + vfTails + vfHeadR
+
+    # All the volume that's left, we will fill with water
+    vfWat = 1 - vfTot
+
+    # Now convert all the Volume Fractions to SLDs
+    sld_Value_Tails = tails / vTail
+    sld_Value_Head = heads / vHead
+
+    sldSilicon = vfSilicon * 2.073e-6
+    sldOxide = vfOxide * 3.41e-6
+
+    sldHeadL = vfHeadL * sld_Value_Head
+    sldHeadR = vfHeadR * sld_Value_Head
+    sldTails = vfTails * sld_Value_Tails
+    sldWat = vfWat * bulk_out[contrast]
+
+    # Put this all together
+    totSLD = sldSilicon + sldOxide + sldHeadL + sldTails + sldHeadR + sldWat
+
+    # Make the SLD array for output
+    SLD = [[a, b] for (a, b) in zip(z, totSLD)]
+
+    return SLD, subRough
+
+
+def layer(z, prevLaySurf, thickness, height, Sigma_L, Sigma_R):
+    """This produces a layer, with a defined thickness, height and roughness.
+    Each side of the layer has its own roughness value.
+    """
+    # Find the edges
+    left = prevLaySurf
+    right = prevLaySurf + thickness
+
+    # Make our heaviside
+    a = (z - left) / ((2**0.5) * Sigma_L)
+    b = (z - right) / ((2**0.5) * Sigma_R)
+
+    erf_a = np.array([math.erf(value) for value in a])
+    erf_b = np.array([math.erf(value) for value in b])
+
+    VF = np.array((height / 2) * (erf_a - erf_b))
+
+    return VF, right

examples/DSPC_custom_layers/custom_bilayer_DSPC.py

@@ -0,0 +1,85 @@
+import numpy as np
+
+
+def custom_bilayer_DSPC(params, bulk_in, bulk_out, contrast):
+    """CUSTOMBILAYER RAT Custom Layer Model File.
+
+    This file accepts 3 vectors containing the values for params, bulk in and bulk out.
+    The final parameter is an index of the contrast being calculated.
+
+    The function should output a matrix of layer values, in the form...
+
+    Output = [thick 1, SLD 1, Rough 1, Percent Hydration 1, Hydrate how 1
+              ....
+              thick n, SLD n, Rough n, Percent Hydration n, Hydration how n]
+
+    The "hydrate how" parameter decides if the layer is hydrated with Bulk out or Bulk in phases.
+    Set to 1 for Bulk out, zero for Bulk in.
+    Alternatively, leave out hydration and just return...
+
+    Output = [thick 1, SLD 1, Rough 1,
+              ....
+              thick n, SLD n, Rough n]
+
+    The second output parameter should be the substrate roughness.
+    """
+    sub_rough = params[0]
+    oxide_thick = params[1]
+    oxide_hydration = params[2]
+    lipidAPM = params[3]
+    headHydration = params[4]
+    bilayerHydration = params[5]
+    bilayerRough = params[6]
+    waterThick = params[7]
+
+    # We have a constant SLD for the bilayer
+    oxide_SLD = 3.41e-6
+
+    # Now make the lipid layers
+    # Use known lipid volume and compositions to make the layers
+
+    # define all the neutron b's.
+    bc = 0.6646e-4  # Carbon
+    bo = 0.5843e-4  # Oxygen
+    bh = -0.3739e-4  # Hydrogen
+    bp = 0.513e-4  # Phosphorus
+    bn = 0.936e-4  # Nitrogen
+
+    # Now make the lipid groups
+    COO = (4 * bo) + (2 * bc)
+    GLYC = (3 * bc) + (5 * bh)
+    CH3 = (2 * bc) + (6 * bh)
+    PO4 = (1 * bp) + (4 * bo)
+    CH2 = (1 * bc) + (2 * bh)
+    CHOL = (5 * bc) + (12 * bh) + (1 * bn)
+
+    # Group these into heads and tails:
+    Head = CHOL + PO4 + GLYC + COO
+    Tails = (34 * CH2) + (2 * CH3)
+
+    # We need volumes for each. Use literature values:
+    vHead = 319
+    vTail = 782
+
+    # We use the volumes to calculate the SLDs
+    SLDhead = Head / vHead
+    SLDtail = Tails / vTail
+
+    # We calculate the layer thickness' from the volumes and the APM
+    headThick = vHead / lipidAPM
+    tailThick = vTail / lipidAPM
+
+    # Manually deal with hydration for layers in this example.
+    oxSLD = (oxide_hydration * bulk_out[contrast]) + ((1 - oxide_hydration) * oxide_SLD)
+    headSLD = (headHydration * bulk_out[contrast]) + ((1 - headHydration) * SLDhead)
+    tailSLD = (bilayerHydration * bulk_out[contrast]) + ((1 - bilayerHydration) * SLDtail)
+
+    # Make the layers
+    oxide = [oxide_thick, oxSLD, sub_rough]
+    water = [waterThick, bulk_out[contrast], bilayerRough]
+    head = [headThick, headSLD, bilayerRough]
+    tail = [tailThick, tailSLD, bilayerRough]
+
+    output = np.array([oxide, water, head, tail, tail, head])
+
+    return output, sub_rough

examples/absorption/volume_thiol_bilayer.py

@@ -0,0 +1,138 @@
+def volume_thiol_bilayer(params, bulk_in, bulk_out, contrast):
+    """VolumeThiolBilayer  RAT Custom Layer Model File.
+
+    This file accepts 3 vectors containing the values for params, bulk in and bulk out.
+    The final parameter is an index of the contrast being calculated
+
+    The function should output a matrix of layer values, in the form...
+
+    Output = [thick 1, SLD 1, Rough 1, Percent Hydration 1, Hydrate how 1
+              ....
+              thick n, SLD n, Rough n, Percent Hydration n, Hydration how n]
+
+    The "hydrate how" parameter decides if the layer is hydrated with Bulk out or Bulk in phases.
+    Set to 1 for Bulk out, zero for Bulk in.
+    Alternatively, leave out hydration and just return...
+
+    Output = [thick 1, SLD 1, Rough 1,
+              ....
+              thick n, SLD n, Rough n]
+
+    The second output parameter should be the substrate roughness.
+    """
+    subRough = params[0]
+    alloyThick = params[1]
+    alloySLDUp = params[2]
+    alloyISLDUp = params[3]
+    alloySLDDown = params[4]
+    alloyISLDDown = params[5]
+    alloyRough = params[6]
+    goldThick = params[7]
+    goldRough = params[8]
+    goldSLD = params[9]
+    goldISLD = params[10]
+    thiolAPM = params[11]
+    thiolHeadHydr = params[12]
+    thiolCoverage = params[13]
+    cwThick = params[14]
+    bilayerAPM = params[15]
+    bilHeadHydr = params[16]
+    bilayerRough = params[17]
+    bilayerCoverage = params[18]
+
+    # Make the metal layers
+    gold = [goldThick, goldSLD, goldISLD, goldRough]
+    alloyUp = [alloyThick, alloySLDUp, alloyISLDUp, alloyRough]
+    alloyDown = [alloyThick, alloySLDDown, alloyISLDDown, alloyRough]
+
+    # Neutron b's..
+    # define all the neutron b's.
+    bc = 0.6646e-4  # Carbon
+    bo = 0.5843e-4  # Oxygen
+    bh = -0.3739e-4  # Hydrogen
+    bp = 0.513e-4  # Phosphorus
+    bn = 0.936e-4  # Nitrogen
+
+    # Work out the total scattering length in each fragment
+    # Define scattering lengths
+    # Hydrogenated version
+    COO = (2 * bo) + (bc)
+    GLYC = (3 * bc) + (5 * bh)
+    CH3 = (1 * bc) + (3 * bh)
+    PO4 = (1 * bp) + (4 * bo)
+    CH2 = (1 * bc) + (2 * bh)
+    CH = (1 * bc) + (1 * bh)
+    CHOL = (5 * bc) + (12 * bh) + (1 * bn)
+
+    # And also volumes
+    vCH3 = 52.7  # CH3 volume in the paper appears to be for 2 * CH3's
+    vCH2 = 28.1
+    vCOO = 39.0
+    vGLYC = 68.8
+    vPO4 = 53.7
+    vCHOL = 120.4
+    vCHCH = 42.14
+
+    vHead = vCHOL + vPO4 + vGLYC + 2 * vCOO
+    vTail = (28 * vCH2) + (1 * vCHCH) + (2 * vCH3)  # Tail volume
+
+    # Calculate sum_b's for other fragments
+    sumbHead = CHOL + PO4 + GLYC + 2 * COO
+    sumbTail = (28 * CH2) + (2 * CH) + 2 * CH3
+
+    # Calculate SLDs and Thickness
+    sldHead = sumbHead / vHead
+    thickHead = vHead / thiolAPM
+
+    sldTail = sumbTail / vTail
+    thickTail = vTail / thiolAPM
+
+    # Correct head SLD based on hydration
+    thiolHeadHydr = thiolHeadHydr / 100
+    sldHead = sldHead * (1 - thiolHeadHydr) + (thiolHeadHydr * bulk_out[contrast])
+
+    # Now correct both the SLDs for the coverage parameter
+    sldTail = (thiolCoverage * sldTail) + ((1 - thiolCoverage) * bulk_out[contrast])
+    sldHead = (thiolCoverage * sldHead) + ((1 - thiolCoverage) * bulk_out[contrast])
+
+    SAMTAILS = [thickTail, sldTail, 0, goldRough]
+    SAMHEAD = [thickHead, sldHead, 0, goldRough]
+
+    # Now do the same for the bilayer
+    vHead = vCHOL + vPO4 + vGLYC + 2 * vCOO
+    vTail = 28 * vCH2  # Tail volume
+    vMe = 2 * vCH3
+
+    sumbHead = CHOL + PO4 + GLYC + 2 * COO
+    sumbTail = 28 * CH2
+    sumbMe = 2 * CH3
+
+    sldHead = sumbHead / vHead
+    thickHead = vHead / bilayerAPM
+    bilHeadHydr = bilHeadHydr / 100
+    sldHead = sldHead * (1 - bilHeadHydr) + (bilHeadHydr * bulk_out[contrast])
+
+    sldTail = sumbTail / vTail
+    thickTail = vTail / bilayerAPM
+
+    sldMe = sumbMe / vMe
+    thickMe = vMe / bilayerAPM
+
+    sldTail = (bilayerCoverage * sldTail) + ((1 - bilayerCoverage) * bulk_out[contrast])
+    sldHead = (bilayerCoverage * sldHead) + ((1 - bilayerCoverage) * bulk_out[contrast])
+    sldMe = (bilayerCoverage * sldMe) + ((1 - bilayerCoverage) * bulk_out[contrast])
+
+    BILTAILS = [thickTail, sldTail, 0, bilayerRough]
+    BILHEAD = [thickHead, sldHead, 0, bilayerRough]
+    BILME = [thickMe, sldMe, 0, bilayerRough]
+
+    BILAYER = [BILHEAD, BILTAILS, BILME, BILME, BILTAILS, BILHEAD]
+
+    CW = [cwThick, bulk_out[contrast], 0, bilayerRough]
+
+    if contrast == 1 or contrast == 3:
+        output = [alloyUp, gold, SAMTAILS, SAMHEAD, CW, *BILAYER]
+    else:
+        output = [alloyDown, gold, SAMTAILS, SAMHEAD, CW, *BILAYER]
+
+    return output, subRough