Creating Mathematica Functions to Determine Degree of Crystallinity from XRD Plots

Building a Manual XRD Curve Fitting Program to Calculate Degree of Crystallinity

Summary

We can calculate percent crystallinity of semi-crystalline polymers from XRD plots. These XRD plots show crystalline peaks and amorphous peaks.  Taking the ratio of the area under the crystalline peaks over the total area under all peaks gives us the degree of crystallinity.  Peak fitting software, or in our case Mathematica,  can be used to identify and perform calculations on these peaks.

Generating a Manipulable Gaussian Peak in Mathematica

The normal probability density function can be used to generate a gaussian peak.  In this function, sigma is standard deviation and mu is mean.

Here I define a function for a gaussian curve in Mathematica.  Sigma is replaced with width and mu is replaced with shift.  An additional variable, height, allows me to change the curve independently of the height defined by a given standard deviation.  Using this code, I can create peaks of varying widths, heights, and x positions.

gaussianCurve=1/(width*Sqrt[2 Pi])*Exp[(-(x-shift)^2)/(2 width^2)]*height;

With[{localGaussianCurve=gaussianCurve},
Manipulate[Plot[localGaussianCurve,{x,-10,10},PlotStyle->Red,PlotRange->{{-10,10},{0,1}}],{width,0.01,5},{shift,-10,10},{height,0,10}]]

Utilizing Multiple Manipulable Peaks to Create a Curve Fitting Program

Most XRD curves will have multiple peaks to fit.  Below I define 3 different peaks.

curveOne=1/(width1*Sqrt[2 Pi])*Exp[(-(x-shift1)^2)/(2 width1^2)]*height1;

curveTwo=1/(width2*Sqrt[2 Pi])*Exp[(-(x-shift2)^2)/(2 width2^2)]*height2;

curveThree=1/(width3*Sqrt[2 Pi])*Exp[(-(x-shift3)^2)/(2 width3^2)]*height3;

Here I created a graph with three peaks.  The user can manipulate peaks and choose which peaks to display.  A orange line shows the combined values of all three peaks.

combinedCurves =With[{localCurveOne=curveOne,localCurveTwo=curveTwo,localCurveThree=curveThree,total=curveOne+curveTwo+curveThree},Manipulate[Plot[{Evaluate[peakNum/.{1->localCurveOne,2->localCurveTwo,3->localCurveThree}],{localCurveOne+localCurveTwo+localCurveThree}},{x,-10,10},PlotRange->{{-10,10},{0,1}}],{width1,0.01,5},{height1,0,5},{shift1,-10,10},{width2,0.01,5},{height2,0,5},{shift2,-10,10},{width3,0.01,5},{height3,0,5},{shift3,-10,10},{{peakNum,{1,2,3}},{1,2,3},CheckboxBar}]]

Here I uploaded a real semi-crystalline XRD plot to match.  The user can manipulate both the size of the XRD plot and the location for better viewing.  I also added in the ability to save the equations of each peak.

mySampleCurve=;

combinedCurves =With[{localCurveOne=curveOne,localCurveTwo=curveTwo,localCurveThree=curveThree},Manipulate[Column[{Plot[{Evaluate[peakNum/.{1->localCurveOne,2->localCurveTwo,3->localCurveThree}],{localCurveOne+localCurveTwo+localCurveThree}},{x,-10,10},PlotRange->{{-10,10},{0,1}},Prolog->{Inset[mySampleCurve,pos,{0,0},scale]},PlotRangeClipping->False,AxesStyle->Black],Row[{"Crystalline Peak 1",Set[crystallinePeak1LDPE,localCurveOne]}],Row[{"Crystalline Peak 2",Set[crystallinePeak2LDPE,localCurveTwo]}],Row[{"Amorphous Peak",Set[amorphousPeak,localCurveThree]}]}],{{pos,{0,0}},{-15,-2},{2,2}},Dynamic@pos,{{scale,10},1,24},Delimiter,Style["Crystalline Peak 1",12,Bold],{{width1,1,"Width"},0.01,5},{{height1,1,"Height"},0,5},{{shift1,1,"Shift"},-10,10},Delimiter,Style["Crystalline Peak 2",12,Bold],{{width2,1,"Width"},0.01,5},{{height2,1,"Height"},0,5},{{shift2,1,"Shift"},-10,10},Delimiter,Style["Amorphous Peak",12,Bold],{{width3,1,"Width"},0.01,5},{{height3,1,"Height"},0,5},{{shift3,1,"Shift"},-10,10},{{peakNum,{1,2,3}},{1,2,3},CheckboxBar},SaveDefinitions->True]]

Calculating Degree of Crystallinity

degreeOfCyrstallinity=(Integrate[crystallinePeak1LDPE,{x,-10,10}]+Integrate[crystallinePeak2LDPE,{x,-10,10}])/(Integrate[crystallinePeak1LDPE,{x,-10,10}]+Integrate[crystallinePeak2LDPE,{x,-10,10}]+Integrate[amorphousPeak,{x,-10,10}])

Polyethylene Case Study

All of the plots below were generated using the Mathematica peak fitting program described above!

The expected values for degree of crystallinity are as follows:
Low Density Polyethylene (LDPE): 44-51% crystallinity 
High Density Polyethylene(HDPE): 48-69% crystallinity 
Ultra High Molecular Weight Polyethylene: 75-77% crystallinity

Low Density Polyethylene

Mathematica Fitting of a Low Density Polyethylene XRD Pattern

Calculating Degree of Crystallinity for Low Density Polyethylene

The expected value for degree of crystallinity for low density polyethylene (LDPE) is 44-51% crystallinity.  The value obtained using Mathematica peak fitting is 46% crystallinity.  This is within the accepted range.

High Density Polyethylene

Mathematica Fitting of a High Density Polyethylene XRD Pattern

Calculating Degree of Crystallinity for High Density Polyethylene

The expected value for degree of crystallinity for high density polyethylene (HDPE) is 48-69% crystallinity.  The value obtained using Mathematica peak fitting is 56% crystallinity.  This is within the accepted range.

Ultra High Molecular Weight Polyethylene

Mathematica Fitting of a Ultra High Molecular Weight Polyethylene XRD Pattern

Calculating Degree of Crystallinity for Ultra High Molecular Weight Polyethylene

The expected value for degree of crystallinity for ultra high molecular weight polyethylene (UHMW) is 75-77% crystallinity.  The value obtained using Mathematica peak fitting is 76% crystallinity.  This is within the accepted range.

Understanding Structural Properties and Applications of Polyethylene

LDPE is branched while HDPE is more linear. As such, HDPE is able to pack more tightly than LDPE and has a higher % crystallinity.  LDPE is more flexible than HDPE.  LDPE is used in applications such as plastic bags, plastic cling wrap, and packaging.  HDPE is found in milk jugs, water bottles, and children’s toys. 

UHMW Polyethylene is made up of extremely long chains of polyethylene, which all align in the same direction. It derives its strength largely from the length of each individual molecule (chain). Van der Waals bonds between individual molecules are relatively weak, but because the chains are very long, large overlaps can exist and the material is able to carry large shear forces.  When formed into fibers, the polymer chains can attain a parallel orientation and a degree of crystallinity up to 77% crystallinity.  Currently UHMW Polyethylene has the highest impact strength of any thermoplastic.

Building an Advanced Curve Fitting Program to Calculate Degree of Crystallinity

Summary

We can calculate percent crystallinity of semi-crystalline polymers from XRD plots. These XRD plots show crystalline peaks and amorphous peaks.  Taking the ratio of the area under the crystalline peaks over the total area under all peaks gives us the degree of crystallinity.  Peak fitting software, or in our case Mathematica,  can be used to identify and perform calculations on these peaks.  Previously, this notebook demonstrated a peak fitting program that used user guesses to estimate degree of crystallinity.  The code described below will allow the user to make a first guess of peak position, and Mathematica will help to find the best fit for the data.

Extracting Data from Images

Here I import the image I wish to extract data from.

LDPEClean=;

My original image was blue and black. This code makes the pixels either white or black.

curve=ImageApply[{0,0,0}&,LDPEClean,Masking->ColorNegate[Binarize[GaussianFilter[LDPEClean,5]]]]

I extract the position of all of the black pixels.

curvLoc=(Reverse/@Position[ImageData[curve,DataReversed->True],{0.,0.,0.}])

ListPlot[curvLoc, PlotRange->Full]

The new data plots a rather thick line.  We want to know a more exact position for the curve, so we can find the average y value of the curve at each x position.

singleCurveLoc=Table[Mean[Cases[curvLoc,{i,_}]],{i,1,Max[curvLoc]}];

ListPlot[singleCurveLoc,PlotRange->All]

The averaged y value predicts peak heights that are slightly too low for narrow peaks due to the thickness of the plot line.  To fix this problem, we can extract either the greatest or second greatest y value.

singleCurveLocMax=Table[{i,Max[#[[2]]&/@Cases[curvLoc,{i,_}]]},{i,1,Max[curvLoc]}];

ListPlot[singleCurveLocMax, PlotRange->All]

singleCurveLocSecondMax=Table[{i,RankedMax[#[[2]]&/@Cases[curvLoc,{i,_}],2]},{i,1,Max[curvLoc]}];

ListPlot[singleCurveLocSecondMax, PlotRange->All]

background=ListPlot[singleCurveLocSecondMax, PlotRange->All,ImageSize->Large]

Utilizing User Guesses for Peaks and Mathematic Best Fit Functions to Create a Curve Fitting Program

Like before, this program will try to match the provided data to three peaks.  The model below allows height of peaks, width of peaks, and peak x position to be manipulated.

model=height+amp1*Exp[-(x-x01)^2/sigma1^2]+amp2*Exp[-(x-x02)^2/sigma2^2]+amp3*Exp[-(x-x03)^2/sigma3^2]

This code allows the user to manipulate all three curves and put in best guesses for peak positions.  Using these guesses, Mathematica will find a set of three gaussian curves that are the closest match.

findBestFitFromValues[{amp1guess_,x01guess_,sigma1guess_,amp2guess_,x02guess_,sigma2guess_,amp3guess_,x03guess_,sigma3guess_,heightguess_}]:=
FindFit[singleCurveLocSecondMax,{model,{sigma1>0,sigma2>0,sigma3>0}},{{amp1,amp1guess},{x01,x01guess},{sigma1,sigma1guess},{amp2,amp2guess},{x02,x02guess},{sigma2,sigma2guess},{amp3,amp3guess},{x03,x03guess},{sigma3,sigma3guess},{height,heightguess}},x](*This is a function that takes guesses and finds the best fit. Sigma was constrained to be positive.*)

With[
{
localModel=
model/.
{
amp1-> amp1Var, amp2-> amp2Var,amp3-> amp3Var,
sigma1->sigma1Var,sigma2-> sigma2Var,sigma3-> sigma3Var,
x01-> x01Var, x02-> x02Var,x03-> x03Var,
height-> heightVar
}},
Manipulate[
Column[{
Style["Match to Data",12,Bold],
Show[background,Plot[localModel,{x,0,500},PlotRange->All],Graphics[
{
Orange,Line[{{x01Var,0}, {x01Var,500}}],
Blue,Line[{{x02Var,0}, {x02Var,500}}],
Red,Line[{{x03Var,0}, {x03Var,500}}]
}
]],
Style["Chosen Curve",12,Bold],
Plot[localModel,{x,0,500},PlotRange->Full]}
],
Delimiter,Style["Crystalline Peak 1",12,Bold],
{{amp1Var,225, Style["Amplitude 1", Orange]},0,400},
{{x01Var,250,Style["center 1", Orange]},0,400},
{{sigma1Var,10,Style["sigma 1", Orange]},0,100},
Delimiter,Style["Crystalline Peak 2",12,Bold],
{{amp2Var,60,Style["Amplitude 2", Blue]},0,400},
{{x02Var,290,Style["center 2", Blue]},0,500},
{{sigma2Var,15,Style["sigma 2", Blue]},0,70},
Delimiter,Style["Amorphous Peak",12,Bold],
{{amp3Var,45,Style["Amplitude 3", Red]},0,400},
{{x03Var,210,Style["center 3", Red]},0,500},
{{sigma3Var,50,Style["sigma 3", Red]},0,100},
Delimiter,Style["Height",12,Bold],
{{heightVar,15,Style["Height"]},0,50},
Delimiter,
Control[Button["click find nearest solution",
vals =
{amp1Var,x01Var,sigma1Var,amp2Var,x02Var,sigma2Var,amp3Var,x03Var,sigma3Var,heightVar} ={amp1,x01,sigma1,amp2,x02,sigma2,amp3,x03,sigma3,height}/.
findBestFitFromValues[
{amp1Var,x01Var,sigma1Var,amp2Var,x02Var,sigma2Var,amp3Var,x03Var,sigma3Var,heightVar}]]],Delimiter,Style["Obtained Values",12,Bold],
Row[{
Dynamic[
{
Set[amp1UserDefined,amp1Var],
Set[x01UserDefined,x01Var],
Set[sigma1UserDefined,sigma1Var],
Set[amp2UserDefined,amp2Var],
Set[x02UserDefined,x02Var],
Set[sigma2UserDefined,sigma2Var],
Set[amp3UserDefined,amp3Var],
Set[x03UserDefined,x03Var],
Set[sigma3UserDefined,sigma3Var],
Set[heightUserDefined,heightVar]},"  "
]}],
SaveDefinitions->True
]
]

Calculating Degree of Crystallinity

curve1=Integrate[amp1UserDefined*Exp[-(x-x01UserDefined)^2/sigma1UserDefined^2],{x,0,450}];
curve2=Integrate[amp2UserDefined*Exp[-(x-x02UserDefined)^2/sigma2UserDefined^2],{x,0,450}];
curve3=Integrate[amp3UserDefined*Exp[-(x-x03UserDefined)^2/sigma3UserDefined^2],{x,0,450}];
DegreeOfCrystallinityNew=(curve1+curve2)/(curve1+curve2+curve3)