Cross-correlated relaxation in the triyne molecule

List of symbols

DD: “dipole-dipole” - e.g. the “DDtensor” is the intrapair dipole-dipole tensor
bjk: intrapair dipolar coupling constant for spins j and k
bCC: bjk, where j and k are the13C nuclei in question
δ: chemical shift - e.g. “δtensor” is the chemical shift tensor
CSA: chemical shift anisotropy
Atensor: irreducible spherical tensor form of the corresponding Cartesian spatial tensor
AAS: arbitrary axis system (that which the computational results are expressed in)
P, P-frame: principal axis system of a given interaction
D, D-frame: diffusion frame
Hcoh: coherent Hamiltonian
Γ: relaxation superoperator
SC: semi-classical
Lb: Lindbladian
ωΔ: Difference in chemically-shifted Larmor frequencies
ωΣ: sum of chemically-shifted Larmor frequencies
ωJ: 2(Subscript[πJ, jk]), where (Subscript[J, jk])is the isotropic spin-spin coupling

Initialisation

Load SpinDynamica

Needs["SpinDynamica`"]

Consider 2 spin-1/2 nuclei

SetSpinSystem[2]

Derive DD and CSA spatial tensors in P- (or D-) frame

Form Cartesian DD tensor in the principal axis system (P-frame)

DDCartesianTensor[{"C1","C2"},"P"]=bjk*{{-1,0,0},{0,-1,0},{0,0,2}};

DDCartesianTensor[{"C1","C2"},"P"]//MatrixForm

Transform Cartesian tensor to irreducible spherical tensor

Atensor[2,{"DD",{"C1","C2"}},"P"]=CartesianToSpherical[2,DDCartesianTensor[{"C1","C2"},"P"]]

Atensor[2,{"DD",{"C1","C2"}},"P"][0]

Form rank 0 and 2 chemical shift Cartesian tensors in the P-frame

δtensor["C1","rank0","P"]=(Subscript[δ, iso])["C1"]IdentityMatrix[3];

δtensor["C1","rank0","P"]//MatrixForm

δtensor["C2","rank0","P"]=(Subscript[δ, iso])["C2"]IdentityMatrix[3];

δtensor["C2","rank0","P"]//MatrixForm

δtensor["C1","rank2","P"]=δaniso["C1"]{{-(1/2)(1-η["C1"]),0,0},{0,-(1/2)(1+η["C1"]),0},{0,0,1}};

δtensor["C1","rank2","P"]//MatrixForm

δtensor["C2","rank2","P"]=δaniso["C2"]{{-(1/2)(1-η["C2"]),0,0},{0,-(1/2)(1+η["C2"]),0},{0,0,1}};

δtensor["C2","rank2","P"]//MatrixForm

Put these tensors in irreducible spherical tensor form

δtensor["C1","P"]=δtensor["C1","rank0","P"]+δtensor["C1","rank2","P"];

Atensor[2,{"CSA","C1"},"P"]=CartesianToSpherical[2,δtensor["C1","P"]]//Simplify

δtensor["C2","P"]=δtensor["C2","rank0","P"]+δtensor["C2","rank2","P"];

Atensor[2,{"CSA","C2"},"P"]=CartesianToSpherical[2,δtensor["C2","P"]]//Simplify

Atensor[2,{"CSA","C2"},"P"][0]

Atensor[2,{"CSA","C2"},"P"][-2]

Form spin tensors in L frame

form DD spin tensor in the lab (L-) frame

Xtensor[2,{"DD",{"C1","C2"}},"L"]=SphericalTensor[2,Table[opT[{1,2},{2,m}],{m,2,-2,-1}]];

form CSA spin-field tensors in the L-frame

Btensor[1]=CartesianToSpherical[1,{Bx,By,Bz}]/.{Bx->0,By->0,Bz->B0};

Ttensor[1,{"CS","C1"}]=SphericalTensor[1,Table[opT[1,{1,m}],{m,1,-1,-1}]];

Ttensor[1,{"CS","C2"}]=SphericalTensor[1,Table[opT[2,{1,m}],{m,1,-1,-1}]];

XtensorCSA[2,site_]:=
SphericalTensor[2,
Table[
Sum[CG[{2,m},{1,m1},{1,m2}]×Btensor[1][m1]×Ttensor[1,{"CS",site}][m2],
{m1,-1,1},{m2,-1,1}],
{m,2,-2,-1}]
]

Xtensor[2,{"CSA","C1"},"L"]=XtensorCSA[2,"C1"]

Xtensor[2,{"CSA","C2"},"L"]=XtensorCSA[2,"C2"]

Define spectral density function

j[ω_,τ_]:=τ/(1+(ω^2)(τ^2))

Write (Γ^total) = (Γ^DD) + (Γ^CSA) + (Γ^DDxCSA)

interaction constants

RelaxationConstant[{"CSA","C1"|"C2"}]=-γ;

RelaxationConstant[{"DD",{"C1","C2"}}]=1;

semi-classical expression

RelaxationSuperoperator[interaction1_,interaction2_,w_,τ_]:=-(1/5)RelaxationConstant[interaction1]×RelaxationConstant[interaction2]×Sum[(((-1))^(p+q))Adjoint[Atensor[2,interaction1,"P"]][p]×Atensor[2,interaction2,"P"][q],{p,-2,2},{q,-2,2}]×Sum[DoubleCommutationSuperoperator[Xtensor[2,interaction1,"L"][m],Adjoint[Xtensor[2,interaction2,"L"][m]]]×j[m×w,τ],{m,-2,2}]

ΓDD=RelaxationSuperoperator[{"DD",{"C1","C2"}},{"DD",{"C1","C2"}},w,t];

ΓCSA=RelaxationSuperoperator[{"CSA","C1"},{"CSA","C1"},w,t]+RelaxationSuperoperator[{"CSA","C2"},{"CSA","C2"},w,t]+RelaxationSuperoperator[{"CSA","C1"},{"CSA","C2"},w,t]+RelaxationSuperoperator[{"CSA","C2"},{"CSA","C1"},w,t];

ΓDDxCSA=RelaxationSuperoperator[{"DD",{"C1","C2"}},{"CSA","C1"},w,t]+RelaxationSuperoperator[{"DD",{"C1","C2"}},{"CSA","C2"},w,t]+RelaxationSuperoperator[{"CSA","C1"},{"DD",{"C1","C2"}},w,t]+RelaxationSuperoperator[{"CSA","C2"},{"DD",{"C1","C2"}},w,t];

Γtotal=ΓDD+ΓCSA+ΓDDxCSA;

Lindbladian expression

LindbladianRelaxation[interaction1_,interaction2_,w_,τ_]:=2(1/5)RelaxationConstant[interaction1]×RelaxationConstant[interaction2]×Sum[(((-1))^(p+q))Adjoint[Atensor[2,interaction1,"P"][p]]×Atensor[2,interaction2,"P"][q],{p,-2,2},{q,-2,2}]×Sum[Lindbladian[Xtensor[2,interaction1,"L"][m],Adjoint[Xtensor[2,interaction2,"L"]][m]]×j[m×w,τ],{m,-2,2}]

LbDD=LindbladianRelaxation[{"DD",{"C1","C2"}},{"DD",{"C1","C2"}},w,t];

LbCSA=LindbladianRelaxation[{"CSA","C1"},{"CSA","C1"},w,t]+LindbladianRelaxation[{"CSA","C2"},{"CSA","C2"},w,t]+LindbladianRelaxation[{"CSA","C1"},{"CSA","C2"},w,t]+LindbladianRelaxation[{"CSA","C2"},{"CSA","C1"},w,t];

LbDDxCSA=LindbladianRelaxation[{"DD",{"C1","C2"}},{"CSA","C1"},w,t]+LindbladianRelaxation[{"DD",{"C1","C2"}},{"CSA","C2"},w,t]+LindbladianRelaxation[{"CSA","C1"},{"DD",{"C1","C2"}},w,t]+LindbladianRelaxation[{"CSA","C2"},{"DD",{"C1","C2"}},w,t];

Lbfull=LbDD+LbCSA+LbDDxCSA;

Check routine using (Γ^DD) and (Subscript[B, ST]) basis

SetOperatorBasis[SingletTripletOperatorBasis[]]

SuperoperatorMatrixRepresentation[RelaxationSuperoperator[{"DD",{"C1","C2"}},{"DD",{"C1","C2"}},w,t]]//MatrixForm//FullSimplify;

SuperoperatorMatrixRepresentation[LbDD]//MatrixForm//FullSimplify;

Get analytical matrix elements

define (Subscript[B, ST'])

SetBasis[SingletTripletBasis[]]

S0ket=Cos[θ/2]Ket[1]-Sin[θ/2]Ket[3];

S0bra=Cos[θ/2]Bra[1]-Sin[θ/2]Bra[3];

T0ket=Cos[θ/2]Ket[3]+Sin[θ/2]Ket[1];

T0bra=Cos[θ/2]Bra[3]+Sin[θ/2]Bra[1];

Tp1ket=Ket[2];

Tp1bra=Bra[2];

Tm1ket=Ket[4];

Tm1bra=Bra[4];

define superkets

S0Tp1superket=S0ket.Tp1bra//Simplify

Tm1S0superket=Tm1ket.S0bra//Simplify

T0Tp1superket=T0ket.Tp1bra//Simplify

Tm1T0superket=Tm1ket.T0bra//Simplify

Check all superkets describe (-1)-QC

Ihat=CommutationSuperoperator[opI["z"]];

LiouvilleBracket[S0Tp1superket,Ihat,S0Tp1superket]

LiouvilleBracket[Tm1S0superket,Ihat,Tm1S0superket]

LiouvilleBracket[T0Tp1superket,Ihat,T0Tp1superket]

LiouvilleBracket[Tm1T0superket,Ihat,Tm1T0superket]

ΓDD and LbDD expressions

field=b/.First@Solve[LarmorFrequency["1H",b]/(2π)==-400*(10^6),b]

θST=ArcTan[((LarmorFrequency["13C",field](0.16×(10^(-6))))/(2π 214.15))]

(Subscript[θ, ST])<<1, set (Subscript[θ, ST])->0 for simplicity.

R2["DD","+","SC"]=-LiouvilleBracket[T0Tp1superket,ΓDD,T0Tp1superket,Normalize->True]/.θ->0//FullSimplify

R2["DD","-","SC"]=-LiouvilleBracket[Tm1T0superket,ΓDD,Tm1T0superket,Normalize->True]/.θ->0//FullSimplify

R2DD["DD","+","Lb"]=-LiouvilleBracket[T0Tp1superket,LbDD,T0Tp1superket,Normalize->True]/.θ->0//FullSimplify

R2["DD","-","Lb"]=-LiouvilleBracket[Tm1T0superket,LbDD,Tm1T0superket,Normalize->True]/.θ->0//FullSimplify

(8+29 (t^2) (w^2)+12 (t^4) (w^4))/(1+5 (t^2) (w^2)+4 (t^4) (w^4))//FullSimplify

ΓCSA and LbCSA expressions

Assume axial symmetry; set η[C1] = η[C2] = 0.

Clear[η,δaniso]

R2["CSA","+","SC"]=-LiouvilleBracket[T0Tp1superket,ΓCSA,T0Tp1superket,Normalize->True]/.{θ->0,η["C1"]->0,η["C2"]->0}//FullSimplify

R2["CSA","-","SC"]=-LiouvilleBracket[Tm1T0superket,ΓCSA,Tm1T0superket,Normalize->True]/.{θ->0,Adjoint[δaniso["C1"]]->δaniso["C1"],Adjoint[δaniso["C2"]]->δaniso["C2"],η["C1"]->0,η["C2"]->0}//FullSimplify

R2["CSA","+","Lb"]=-LiouvilleBracket[T0Tp1superket,LbCSA,T0Tp1superket,Normalize->True]/.{θ->0,Adjoint[δaniso["C1"]]->δaniso["C1"],Adjoint[δaniso["C2"]]->δaniso["C2"],η["C1"]->0,η["C2"]->0}//FullSimplify

R2["CSA","-","Lb"]=-LiouvilleBracket[Tm1T0superket,LbCSA,Tm1T0superket,Normalize->True]/.{θ->0,Adjoint[δaniso["C1"]]->δaniso["C1"],Adjoint[δaniso["C2"]]->δaniso["C2"],η["C1"]->0,η["C2"]->0}//FullSimplify

Cross-correlated expressions

R2["CSAxDD","+","SC"]=-LiouvilleBracket[T0Tp1superket,ΓDDxCSA,T0Tp1superket,Normalize->True]/.{θ->0,Adjoint[δaniso["C1"]]->δaniso["C1"],Adjoint[δaniso["C2"]]->δaniso["C2"],η["C1"]->0,η["C2"]->0}//FullSimplify

-((3 B0 bjk t (3+2 (t^2) (w^2)) γ (δaniso["C1"]+δaniso["C2"]))/(20 (1+(t^2) (w^2))))

R2CSAxDDminusSC=-LiouvilleBracket[Tm1T0superket,ΓDDxCSA,Tm1T0superket,Normalize->True]/.{θ->0,Adjoint[δaniso["C1"]]->δaniso["C1"],Adjoint[δaniso["C2"]]->δaniso["C2"],η["C1"]->0,η["C2"]->0}//FullSimplify

R2CSAxDDplusLb=-LiouvilleBracket[T0Tp1superket,LbDDxCSA,T0Tp1superket,Normalize->True]/.{θ->0,Adjoint[δaniso["C1"]]->δaniso["C1"],Adjoint[δaniso["C2"]]->δaniso["C2"]}//FullSimplify

R2CSAxDDminusLb=-LiouvilleBracket[Tm1T0superket,LbDDxCSA,Tm1T0superket,Normalize->True]/.{θ->0,Adjoint[δaniso["C1"]]->δaniso["C1"],Adjoint[δaniso["C2"]]->δaniso["C2"]}//FullSimplify

Derive relaxation parameters from computation results

Use coordinates to obtain the dipole-dipole coupling constant in rad (s^(-1)) and Hz

Coordinates["C1","AAS"]={1.718640,0.028594,0.000311}*(10^(-10));

Coordinates["C2","AAS"]={0.497023,0.088309,0.000345}*(10^(-10));

bCC=DirectDipolarCoupling[13,EuclideanDistance[Coordinates["C1","AAS"],Coordinates["C2","AAS"]]]

bCC/(2π)

Derive CSA and biaxiality from computation of the CS tensor

σtensor["C1","AAS"]={{258.1514,-11.5654,0.0240},
{-10.8106,42.1711,0.0014},{-0.0078,0.0033,38.6536}}*(10^(-6));
MatrixForm[σtensor["C1","AAS"]]

σtensor["C2","AAS"]={{258.0906,-10.9137,-0.0180},
{-10.8451,42.8023,0.0004},{-0.0030,0.0000,38.8968}}*(10^(-6));
MatrixForm[σtensor["C2","AAS"]]

σTMSiso=Mean[{183.9736,183.9809,183.9763,183.9764}*(10^(-6))];

δtensor[site_,frame_]:=-(σtensor[site,frame]-DiagonalMatrix[Table[σTMSiso,{3}]])

δtensor["C1","AAS"];

δtensor["C2","AAS"];

Put the tensors in the P-frame first, the elements will be unordered.

δtensorUnordered["C1","P"]=Inverse[Transpose[Eigenvectors[δtensor["C1","AAS"]]]].δtensor["C1","AAS"].Transpose[Eigenvectors[δtensor["C1","AAS"]]]//Chop;

δtensorUnordered["C2","P"]=Inverse[Transpose[Eigenvectors[δtensor["C2","AAS"]]]].δtensor["C2","AAS"].Transpose[Eigenvectors[δtensor["C2","AAS"]]]//Chop;

MatrixForm[δtensorUnordered["C1","P"]]
MatrixForm[δtensorUnordered["C2","P"]]

Order the elements based on the Haeberlen convention, where |(Subscript[δ, zz])-(Subscript[δ, iso])|≥|(Subscript[δ, xx])-(Subscript[δ, iso])|≥|(Subscript[δ, yy])-(Subscript[δ, iso])|.

δiso["C1"]=(1/3)Tr[δtensor["C1","AAS"]];

Abs[δtensorUnordered["C1","P"][[1,1]]-δiso["C1"]]

Abs[δtensorUnordered["C1","P"][[2,2]]-δiso["C1"]]

Abs[δtensorUnordered["C1","P"][[3,3]]-δiso["C1"]]

δtensor["C1","P"]={{δtensorUnordered["C1","P"][[1,1]],0,0},{0,δtensorUnordered["C1","P"][[2,2]],0},{0,0,δtensorUnordered["C1","P"][[3,3]]}};

δiso["C2"]=(1/3)Tr[δtensor["C2","AAS"]]

Abs[δtensorUnordered["C2","P"][[1,1]]-δiso["C2"]]

Abs[δtensorUnordered["C2","P"][[2,2]]-δiso["C2"]]

Abs[δtensorUnordered["C2","P"][[3,3]]-δiso["C2"]]

δtensor["C2","P"]={{δtensorUnordered["C2","P"][[1,1]],0,0},{0,δtensorUnordered["C2","P"][[2,2]],0},{0,0,δtensorUnordered["C2","P"][[3,3]]}};

δtensor["C1","P"]//MatrixForm
δtensor["C2","P"]//MatrixForm

The biaxiality parameters and CSA are,

η["C1"]=((δtensor["C1","P"][[2,2]]-δtensor["C1","P"][[1,1]])/(δtensor["C1","P"][[3,3]]-δiso["C1"]))

η["C2"]=((δtensor["C2","P"][[2,2]]-δtensor["C2","P"][[1,1]])/(δtensor["C2","P"][[3,3]]-δiso["C2"]))

δaniso["C1"]=δtensor["C1","P"][[3,3]]-δiso["C1"]

δaniso["C2"]=δtensor["C2","P"][[3,3]]-δiso["C2"]

(Subscript[R, 1]) Expression

Get approximate expression for (Subscript[R, 1])

Get R1, which is the inverse of T1.

R1=-LiouvilleBracket[opI[1,"z"],Γtotal,opI[1,"z"],Normalize->True]/.{w->LarmorFrequency["13C",field],γ->MagnetogyricRatio["13C"],bjk->bCC,B0->field}//FullSimplify

T1=1/R1

Solve for (Subscript[τ, ⊥]) to estimate correlation time using the experimental (Subscript[T, 1]).

τperp=t/.First@NSolve[T1==2.235,t,Reals]

Coherent Hamiltonian

Get eigenfrequencies of [ Hcoh , ... ]

Hcoh=(1/2)ωΣ(opI[1,"z"]+opI[2,"z"])+(1/2)ωΔ(opI[1,"z"]-opI[2,"z"])+ωJ opI[1].opI[2];

HcohSuperOp=CommutationSuperoperator[Hcoh]

-LiouvilleBracket[S0Tp1superket,HcohSuperOp,S0Tp1superket]/.θ->ArcTan[ωΔ/(ωJ)]//FullSimplify

-LiouvilleBracket[Tm1S0superket,HcohSuperOp,Tm1S0superket]/.θ->ArcTan[ωΔ/(ωJ)]//FullSimplify

-LiouvilleBracket[T0Tp1superket,HcohSuperOp,T0Tp1superket]/.θ->ArcTan[ωΔ/(ωJ)]//FullSimplify

%/.{ωJ->2π*214,ωΔ->LarmorFrequency["13C",field]*(0.16*(10^(-6))),ωΣ->0}

-LiouvilleBracket[Tm1T0superket,HcohSuperOp,Tm1T0superket]/.θ->ArcTan[ωΔ/(ωJ)]//FullSimplify

%/.{ωJ->2π*214,ωΔ->LarmorFrequency["13C",field]*(0.16*(10^(-6))),ωΣ->0}

Plot spectrum

Define the coherent Hamiltonian in the rotating frame with empirical parameters

Hcoh=(1/2)ωΣ(opI[1,"z"]+opI[2,"z"])+(1/2)ωΔ(opI[1,"z"]-opI[2,"z"])+ωJ opI[1].opI[2]/.{ωΣ->0,ωΔ->LarmorFrequency["13C",field](0.16×(10^(-6))),ωJ->2π 214.15}

Define some relaxation parameters

parameters={bjk->bCC,w->LarmorFrequency["13C",field],t->τperp,γ->MagnetogyricRatio["13C"],B0->field}

spectrum=ListPlot[Re@FT@Signal1D[{{2π 1000,"16k"}},BackgroundGenerator->CombineGenerators[Hcoh,Γtotal/.parameters],LineBroadening->2π*0.15],PlotRange->{{-5,5},All},PlotStyle->{Black,Thick},Frame->{False,False,False,False}]

(Subscript[R, 2]) against field plots

Parameters

ParametersAnalytical={θ->0,bjk->-4152.84×2π,t->τperp,δ->-145.7*(10^(-6)),γ->MagnetogyricRatio["13C"],τ->τperp};

Define (Subscript[R, 2]) components symbolically

Clear[δaniso,η]

R2["+"]=-LiouvilleBracket[T0Tp1superket,Γtotal,T0Tp1superket,Normalize->True]/.{θ->0,η["C1"]->0,η["C2"]->0,Adjoint[δaniso["C1"]]->δaniso["C1"],Adjoint[δaniso["C2"]]->δaniso["C2"],w->(-γ B0),t->τ}//FullSimplify

Behaviour in the extreme-narrowing limit

ExtremeNarrow["+","Symbolic"]=(6/5)(bjk^2)τ+(9/10)bjk(-GyromagneticRatio["13C"]*B0)δ τperp+(13/20)(((-GyromagneticRatio["13C"]*B0))^2)(δ^2)τ;

ExtremeNarrow["+"]=ExtremeNarrow["+","Symbolic"]/.ParametersAnalytical

ExtremeNarrow["-","Symbolic"]=(6/5)(bjk^2)τ-(9/10)bjk(-γ B0)δ τperp+(13/20)(((-γ B0))^2)(δ^2)τ

ExtremeNarrow["-"]=ExtremeNarrow["-","Symbolic"]/.ParametersAnalytical

FindMinimum[ExtremeNarrow["+"],B0]

Plot[{ExtremeNarrow["+"],ExtremeNarrow["-"]},{B0,0,8},PlotStyle->{{Black,Thick},{Black,Thick,Dashed}},FrameStyle->Directive[Black,Thick],Frame->{True,True,False,False},LabelStyle->Directive[Black,FontFamily->"Times",FontSize->40],PlotRange->{{0,8},{0,0.5}},PlotLegends->Placed[{(SuperPlus[(Subscript[R, 2])]),(SuperMinus[(Subscript[R, 2])])},{0.25,0.75}],FrameTicks->{Automatic,{0.1,0.2,0.3,0.4,0.5}}]

Behaviour in the long-(Subscript[τ, ⊥]) limit

longτ["-"]=(1/20)(((3 bjk-2 (-γ B0)δ))^2)/.ParametersAnalytical

FindMinimum[longτ["-"],B0]

longτ["+"]=(1/20)(((3 bjk+2 (-γ B0)δ))^2)/.ParametersAnalytical

Plot[{longτ["+"]/(10^9),longτ["-"]/(10^9)},{B0,0,8},PlotStyle->{{Black,Thick},{Black,Thick,Dashed}},FrameStyle->Directive[Black,Thick],Frame->{True,True,False,False},LabelStyle->Directive[Black,FontFamily->"Times",FontSize->40],PlotRange->{{0,8},{0,2.5}},FrameTicks->{Automatic,{0.5,{1.0,"1.0"},1.5,{2.0,"2.0"},2.5}},PlotLegends->Placed[{(SuperPlus[(Subscript[R, 2])]),(SuperMinus[(Subscript[R, 2])])},{0.25,0.75}]]