/******************************************************************************* * McStas, neutron ray-tracing package * Copyright 1997-2003, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: SANSNanodiscsWithTags * * %I * Written by: Martin Cramer Pedersen (mcpe@nbi.dk) * Date: October 17, 2012 * Origin: KU-Science * * A sample of monodisperse phospholipid bilayer nanodiscs in solution (water) - with * histidine tag still on the belt proteins. * * %D * A component simulating the scattering from a box-shaped, thin solution (water) * of monodisperse phospholipid bilayer nanodiscs - with histidine tags still * on the belt proteins. * * %P * AxisRatio: [] Axis ratio of the bilayer patch. * NumberOfLipids: [] Number of lipids per nanodisc. * AreaPerLipidHeadgroup: [AA^2] Area per lipid headgroup - default is POPC. * HeightOfMSP: [AA] Height of the belt protein - default is MSP1D1. * RadiusOfGyrationForHisTag: [AA] Radius of gyration for the his-tag. * VolumeOfOneMSP: [AA^3] Volume of one belt protein - default is MSP1D1. * VolumeOfHeadgroup: [AA^3] Volume of one lipid headgroup - default is POPC. * VolumeOfCH2Tail: [AA^3] Volume of the CH2-chains of one lipid - default is POPC. * VolumeOfCH3Tail: [AA^3] Volume of the CH3-tails of one lipid - default is POPC. * VolumeOfOneHisTag: [AA^3] Volume of one his-tag. * ScatteringLengthOfOneMSP: [cm] Scattering length of one belt protein - default is MSP1D1. * ScatteringLengthOfHeadgroup: [cm] Scattering length of one lipid headgroup - default is POPC. * ScatteringLengthOfCH2Tail: [cm] Scattering length of the CH2-chains of one lipid - default is POPC. * ScatteringLengthOfCH3Tail: [cm] Scattering length of the CH3-tails of one lipid - default is POPC. * ScatteringLengthOfOneHisTag: [cm] Scattering length of one his-tag. * Roughness: [] Factor used to smear the interfaces of the nanodisc. * Concentration: [mM] Concentration of sample. * RhoSolvent: [cm/AA^3] Scattering length density of solvent - default is D2O. * AbsorptionCrosssection: [1/m] Absorption cross section of the sample. * xwidth: [m] Dimension of component in the x-direction. * yheight: [m] Dimension of component in the y-direction. * zdepth: [m] Dimension of component in the z-direction. * SampleToDetectorDistance: [m] Distance from sample to detector (for focusing the scattered neutrons). * DetectorRadius: [m] Radius of the detector (for focusing the scattered neutrons). * * %E *******************************************************************************/ DEFINE COMPONENT SANSNanodiscsWithTags DEFINITION PARAMETERS () SETTING PARAMETERS (AxisRatio = 1.4, NumberOfLipids = 130.0, AreaPerLipidHeadgroup = 65.0, HeightOfMSP = 24.0, RadiusOfGyrationForHisTag = 10.0, VolumeOfOneMSP = 26296.5, VolumeOfHeadgroup = 319.0, VolumeOfCH2Tail = 818.8, VolumeOfCH3Tail = 108.6, VolumeOfOneHisTag = 2987.3, ScatteringLengthOfOneMSP = 8.8E-10, ScatteringLengthOfHeadgroup = 7.05E-12, ScatteringLengthOfCH2Tail = -1.76E-12, ScatteringLengthOfCH3Tail = -9.15E-13,ScatteringLengthOfOneHisTag = 1.10E-10, Roughness = 5.0, Concentration = 0.01, RhoSolvent = 6.4e-14, AbsorptionCrosssection = 0.0, xwidth, yheight, zdepth, SampleToDetectorDistance, DetectorRadius) OUTPUT PARAMETERS () SHARE %{ // Functions used for compution the formfactor of a cylinder double FormfactorCylinder(double q, double MajorSemiAxis, double MinorSemiAxis, double Height, double Alpha, double Beta) { double const ProjectedRadius = sqrt(pow(MajorSemiAxis * sin(Beta), 2) + pow(MinorSemiAxis * cos(Beta), 2)); double const Formfactor1 = j1(q * ProjectedRadius * sin(Alpha)) / (q * ProjectedRadius * sin(Alpha)); double const Formfactor2 = sin(q * Height * cos(Alpha) / 2.0) / (q * Height * cos(Alpha) / 2.0); return 2 * Formfactor1 * Formfactor2; } // Functions used to compute the characteristics of the histidine tags double Debye(double q, double r) { // Variables needed in functions double ReturnValue; // Constants needed in function double Dummy1; double Dummy2; // Computation Dummy1 = (q * q) * (r * r); Dummy2 = 2 * ((exp(-Dummy1) + Dummy1 - 1) / (Dummy1 * Dummy1)); ReturnValue = Dummy2; return ReturnValue; } double XiRodWithoutEndcaps(double q, double r, double l, double Alpha) { // Variables needed in function double ReturnValue; // Computation ReturnValue = j0(q * r * sin(Alpha)) * sin(q * l * cos(Alpha) / 2.f) / (q * l * cos(Alpha) / 2.f); return ReturnValue; } double PsiHammouda(double x) { // Variables used in function double ReturnValue; // Computation ReturnValue = (1.0 - exp(-pow(x, 2))) / pow(x, 2); return ReturnValue; } // Computation of the intensity from a given nanodisc double IntensityOfEmptyNanodiscsWithTags( double q, double MajorSemiAxis, double MinorSemiAxis, double ThicknessOfBelt, double HeightOfBelt, double HeightOfLipids, double HeightOfTails, double HeightOfCH3, double RadiusOfGyrationForHisTag, double VolumeOfOneHisTag, double DeltaRhoBelt, double DeltaRhoHead, double DeltaRhoCH2Tail, double DeltaRhoCH3Tail, double DeltaRhoTag) { // Declarations double Intensity; double IntensityPart; double IntensityDisc; double AmplitudeOfBelt; double AmplitudeOfHeads; double AmplitudeOfCH2Tail; double AmplitudeOfCH3Tail; double FormfactorOfBelt; double FormfactorOfHeads; double FormfactorOfCH2Tail; double FormfactorOfCH3Tail; const double OuterMajorSemiAxis = MajorSemiAxis + ThicknessOfBelt; const double OuterMinorSemiAxis = MinorSemiAxis + ThicknessOfBelt; const double AverageRadiusOfBelt = sqrt(pow(OuterMajorSemiAxis, 2) + pow(OuterMinorSemiAxis, 2)); const double VolumeOfBelt = PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis; const double VolumeOfHeads = PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis; const double VolumeOfCH2Tail = PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; const double VolumeOfCH3Tail = PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; double Autocorrelation; double Crosscorrelation; double Disccorrelation; double FormfactorOfTags; // Variables needed for integration over alpha int i; const int NumberOfStepsInAlpha = 50; double Alpha; const double AlphaMin = 0.0; const double AlphaMax = PI / 2.0; const double AlphaStep = (AlphaMax - AlphaMin) / (1.0 * NumberOfStepsInAlpha); // Variables needed in integration over beta int j; const int NumberOfStepsInBeta = 50; double Beta; const double BetaMin = 0.0; const double BetaMax = PI / 2.0; const double BetaStep = (BetaMax - BetaMin) / (1.0 * NumberOfStepsInBeta); // Computation Intensity = 0.0; for (i = 0; i < NumberOfStepsInAlpha; ++i) { Alpha = (i + 0.5) * AlphaStep; for (j = 0; j < NumberOfStepsInBeta; ++j) { Beta = (j + 0.5) * BetaStep; // Compute formfactors FormfactorOfBelt = (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis * FormfactorCylinder(q, OuterMajorSemiAxis, OuterMinorSemiAxis, HeightOfBelt, Alpha, Beta) - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis * FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfBelt, Alpha, Beta)) / (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis); FormfactorOfHeads = (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis * FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfLipids, Alpha, Beta) - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails , Alpha, Beta)) / (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis); FormfactorOfCH2Tail = (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta) - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis * FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3 , Alpha, Beta)) / (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis); FormfactorOfCH3Tail = FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta); // Compute amplitudes AmplitudeOfBelt = DeltaRhoBelt * VolumeOfBelt * FormfactorOfBelt; AmplitudeOfHeads = DeltaRhoHead * VolumeOfHeads * FormfactorOfHeads; AmplitudeOfCH2Tail = DeltaRhoCH2Tail * VolumeOfCH2Tail * FormfactorOfCH2Tail; AmplitudeOfCH3Tail = DeltaRhoCH3Tail * VolumeOfCH3Tail * FormfactorOfCH3Tail; // Perform integration IntensityDisc = pow(AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2); // Add the histidine tags FormfactorOfTags = XiRodWithoutEndcaps(q, AverageRadiusOfBelt + RadiusOfGyrationForHisTag, HeightOfBelt, Alpha); Crosscorrelation = 2.0f * pow(DeltaRhoTag * VolumeOfOneHisTag, 2) * pow(FormfactorOfTags, 2) * pow(PsiHammouda(q * RadiusOfGyrationForHisTag), 2); Autocorrelation = 2.0f * pow(DeltaRhoTag * VolumeOfOneHisTag, 2) * Debye(q, RadiusOfGyrationForHisTag); Disccorrelation = 4.0f * (AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail) * DeltaRhoTag * VolumeOfOneHisTag * FormfactorOfTags * PsiHammouda(q * RadiusOfGyrationForHisTag); IntensityPart = IntensityDisc + Crosscorrelation + Autocorrelation + Disccorrelation; // Return the final intensity Intensity += 2 / PI * sin(Alpha) * IntensityPart * AlphaStep * BetaStep; } } return Intensity; } %} DECLARE %{ // Declarations double Absorption; double q; double NumberDensity; // Scattering lengths double RhoBelt; double RhoHead; double RhoCH2Tail; double RhoCH3Tail; double RhoTag; double DeltaRhoHead; double DeltaRhoBelt; double DeltaRhoCH2Tail; double DeltaRhoCH3Tail; double DeltaRhoTag; // Geometric properties double MajorSemiAxis; double MinorSemiAxis; double ThicknessOfBelt; double HeightOfBelt; double HeightOfLipids; double HeightOfTails; double HeightOfCH3; %} INITIALIZE %{ // Rescale concentration into number of aggregates per m^3 times 10^-4 NumberDensity = Concentration * 6.02214129e19; // Computations if (!xwidth || !yheight || !zdepth) { printf("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); } Absorption = AbsorptionCrosssection; // Scattering properties of different components RhoBelt = ScatteringLengthOfOneMSP / VolumeOfOneMSP; RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; RhoTag = ScatteringLengthOfOneHisTag / VolumeOfOneHisTag; DeltaRhoBelt = RhoBelt - RhoSolvent; DeltaRhoHead = RhoHead - RhoSolvent; DeltaRhoCH2Tail = RhoCH2Tail - RhoSolvent; DeltaRhoCH3Tail = RhoCH3Tail - RhoSolvent; DeltaRhoTag = RhoTag - RhoSolvent; // Geometric properties of different components const double AreaOfLipids = NumberOfLipids * AreaPerLipidHeadgroup / 2.0f; MinorSemiAxis = sqrt(AreaOfLipids / (PI * AxisRatio)); MajorSemiAxis = MinorSemiAxis * AxisRatio; HeightOfLipids = 2.0f * (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; HeightOfTails = 2.0f * (VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; HeightOfCH3 = 2.0f * VolumeOfCH3Tail / AreaPerLipidHeadgroup; ThicknessOfBelt = sqrt(pow(MinorSemiAxis + MajorSemiAxis, 2) / 4.0 + 2 * VolumeOfOneMSP / (PI * HeightOfMSP)) - (MajorSemiAxis + MinorSemiAxis) / 2.0; %} TRACE %{ // Declarations double t0; double t1; double l_full; double l; double l1; double Intensity; double Weight; double IntensityPart; double SolidAngle; double qx; double qy; double qz; double v; double dt; double vx_i; double vy_i; double vz_i; char Intersect = 0; // Computation Intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); if (Intersect) { if (t0 < 0.0) { fprintf(stderr, "Neutron already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); ABSORB; } // Compute properties of neutron v = sqrt(pow(vx, 2) + pow(vy, 2) + pow(vz, 2)); l_full = v * (t1 - t0); dt = rand01() * (t1 - t0) + t0; PROP_DT(dt); l = v * (dt - t0); // Store properties of incoming neutron vx_i = vx; vy_i = vy; vz_i = vz; // Generate new direction of neutron randvec_target_circle(&vx, &vy, &vz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); NORM(vx, vy, vz); vx *= v; vy *= v; vz *= v; // Compute q qx = V2K * (vx_i - vx); qy = V2K * (vy_i - vy); qz = V2K * (vz_i - vz); q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); // Compute scattering l1 = v * t1; Intensity = exp(- pow(q * Roughness, 2)) * IntensityOfEmptyNanodiscsWithTags(q, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids, HeightOfTails, HeightOfCH3, RadiusOfGyrationForHisTag, VolumeOfOneHisTag, DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail, DeltaRhoTag); p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp(- Absorption * (l + l1) / v); SCATTER; } %} MCDISPLAY %{ box(0, 0, 0, xwidth, yheight, zdepth); %} END