I'm working on a project for school where I have to develop a large matrix representation of aggregated molecules by focusing on dipole-dipole interactions. Luckily, the aggregate is organized and it produces a symmetric square matrix. I use a sparse array to achieve this. The goal of this project is to predict (and compare to [url=http://dx.doi.org/10.1016/S0006-3495(98)77916-0]these findings) the absorbance and scattering cross-sections of the aggregated molecules when the absorbance of a single molecule is known.
This, in essence, predicts the new absorbance wavelength for the aggregate. Furthermore, given the eigensystem, and our tweaking of Hückel Molecular Orbital Theory that we're using here, we can predict the wavelengths needed to excite electrons in the aggregate to any energy state. If anyone has any questions on the gritty details of the project, please feel free to send me a private message (contact info in my profile).
The simplest system (for a line of molecules) is given below:
n = 36;
sL = SparseArray[{
{i_, i_} -> \[Nu]monL,
{i_, j_} /; Abs[i - j] == 1 -> \[Beta]L},
{n, n}, 0];
My system of functions:
(* Interaction Energies *)
\[Beta][\[Mu]_, r_, \[Theta]_] := (\[Mu]^2*(1 - 3*Cos[\[Theta]]^2))/(4*Pi*c*\[CurlyEpsilon]*h*\[Eta]^2*r^3)
\[Gamma][\[Mu]_, r_, \[Theta]_] := (\[Mu]^2*(1 - 3*Cos[\[Theta]]^2))/(4*Pi*c*\[CurlyEpsilon]*h*\[Eta]^2*r^3)
(* Matrix Operations *)
\[Nu]L[matrix_, minVals_] :=
Eigenvalues[matrix, -minVals, Method -> {"Arnoldi", BasisSize -> n, MaxIterations -> n + 1000}]
\[Mu]L[matrix_, \[Mu]mon_, minVals_] := \[Mu]mon*Abs[Total[Eigenvectors[matrix, -minVals, Method -> {"Arnoldi", BasisSize -> n, MaxIterations -> n + 1000}], {2}]]
(* Energy state values *)
\[Nu]State[eLevel_] := \[Nu]Mins[[Count[\[Nu]Mins, _Real] - (eLevel - 1)]]
\[Mu]State[eLevel_] := \[Mu]Mins[[Count[\[Mu]Mins, _Real] - (eLevel - 1)]]
(* Polarizability functions *)
Re\[Alpha][\[Nu]_, \[Nu]mon_, \[Mu]mon_, eLevel_] := (1/(2*Pi*\[CurlyEpsilon]*h*c))*((\[Nu]mon (\[Nu]mon^2 - \[Nu]^2)*\[Mu]mon^2)/((\[Nu]mon^2 - \[Nu]^2)^2 + \[Nu]^2*\[CapitalGamma]1^2) + Sum[(\[Nu]State[k] (\[Nu]State[k]^2 - \[Nu]^2)*\[Mu]State[k]^2)/((\[Nu]State[k]^2 - \[Nu]^2)^2 + \[Nu]^2*\[CapitalGamma]1^2), {k, 1, eLevel}])
Im\[Alpha][\[Nu]_, \[Nu]mon_, \[Mu]mon_, eLevel_] := (1/(2*Pi*\[CurlyEpsilon]*h*c))*((\[Nu]*\[Nu]mon*\[CapitalGamma]1*\[Mu]mon^2)/((\[Nu]mon^2 - \[Nu]^2)^2 + \[Nu]^2*\[CapitalGamma]1^2) + Sum[(\[Nu]*\[Nu]State[k]*\[CapitalGamma]1*\[Mu]State[k]^2)/((\[Nu]State[k]^2 - \[Nu]^2)^2 + \[Nu]^2*\[CapitalGamma]1^2), {k, 1, eLevel}])
(* Absorbance and Scattering cross-sections *)
Ca[\[Nu]_, \[Nu]mon_, \[Mu]mon_, eLevel_] := 2*Pi*\[Nu]*Im\[Alpha][\[Nu], \[Nu]mon, \[Mu]mon, eLevel]
Cs[\[Nu]_, \[Nu]mon_, \[Mu]mon_, eLevel_] := (8/3)*Pi^3*\[Nu]^4*(Re\[Alpha][\[Nu], \[Nu]mon, \[Mu]mon, eLevel]^2 + Im\[Alpha][\[Nu], \[Nu]mon, \[Mu]mon, eLevel]^2)
Define some constants and variables...
(* Constants *)
\[CurlyEpsilon] = 8.8542*10^(-12);
h = 6.6261*10^(-34);
c = 2.998*10^8;
\[Eta] = 1.33;
\[CapitalGamma]1 = 10000;
\[Lambda]min = 200;
\[Lambda]max = 800;
\[Nu]min = 1/\[Lambda]min*10^9;
\[Nu]max = 1/\[Lambda]max*10^9;
(* Interaction properties *)
r\[Beta] = 0.5*10^(-9);
\[Theta]\[Beta]L = Pi/6;
\[Lambda]L = 435;
\[Mu]monLset = 4;
\[Mu]monL = \[Mu]monLset*3.33564*10^(-30);
\[Nu]monL = 1/\[Lambda]L*10^9;
Compute the interaction energies, and find the 5 smallest eigenvalues and their corresponding eigenvector coefficient totals...
\[Beta]L = \[Beta][\[Mu]monL, r\[Beta], \[Theta]\[Beta]L]
\[Nu]Mins = \[Nu]L[sL, 5]
\[Mu]Mins = \[Mu]L[sL,\[Mu]monL, 5]
And finally, we're able to find some new excitation wavelengths! The lords of physical chemistry be praised!
1/(\[Nu]/.Last[FindMaximum[Ca[\[Nu], \[Nu]monL, \[Mu]monL, 1], {\[Nu], \[Nu]min, \[Nu]max}]])*10^9
1/(\[Nu]/.Last[FindMaximum[Ca[\[Nu], \[Nu]monL, \[Mu]monL, 2], {\[Nu], \[Nu]min, \[Nu]max}]])*10^9
1/(\[Nu]/.Last[FindMaximum[Cs[\[Nu], \[Nu]monL, \[Mu]monL, 1], {\[Nu], \[Nu]min, \[Nu]max}]])*10^9
1/(\[Nu]/.Last[FindMaximum[Cs[\[Nu], \[Nu]monL, \[Mu]monL, 2], {\[Nu], \[Nu]min, \[Nu]max}]])*10^9
The IssueUnfortunately, these maxima are all the same when Mathematica computes them. The results for this particular set-up were computed using MATHCAD around 10 years ago, but without the sparse array ([url=http://dx.doi.org/10.1016/S0006-3495(98)77916-0]the article).
When I try to strong-arm the maximum calculation using Solve[], I get the following:
Solve[D[Ca[\[Nu], \[Nu]monL, \[Mu]monL, 1], \[Nu]] == 0, Reals]
Solve[D[Ca[\[Nu], \[Nu]monL, \[Mu]monL, 2], \[Nu]] == 0, Reals]
Solve[D[Cs[\[Nu], \[Nu]monL, \[Mu]monL, 1], \[Nu]] == 0, Reals]
Solve[D[Cs[\[Nu], \[Nu]monL, \[Mu]monL, 2], \[Nu]] == 0, Reals]
(* Out *)
Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result.
I feel like the issue may lie with the Arnoldi methods that I've set up, but I don't know how to further modify those. I have also tried increasing the working precision in both the eigensystem computation and the maximum computation, but to no avail.
Any help on what I can do to get exact results would be very helpful.
Thank you!
Posting the more instense model in this post would have been a horrifying sight, and would likely have resulted in several people falling out their chairs or throwing their phones and running away from the post. In the article, they performed this analysis for a column-like model. I've constructed a sparse array for a 6 monomer-per-ring model and attached the code.
Edit: bad link
Attachments:
