TensorView.nb

Version 1.2

A simple tool for visualizing second-rank Cartesian tensors within the molecular framework.

Leonard J. Mueller
UC Riverside
April 2013

This tool uses active rotation matrices to manipulate the Cartesian tensor as described in
L.J. Mueller, "Tensors and Rotations in NMR," Concepts in Magnetic Resonance Part A, 38A(5), 221-235 (2011)

Define rotation matrix and fitting function:

Define active rotation matrix

rx[b_]:=(({{1, 0, 0}, {0, Cos[b], -Sin[b]}, {0, Sin[b], Cos[b]}}))
ry[b_]:=(({{Cos[b], 0, Sin[b]}, {0, 1, 0}, {-Sin[b], 0, Cos[b]}}))
rz[b_]:=(({{Cos[b], -Sin[b], 0}, {Sin[b], Cos[b], 0}, {0, 0, 1}}))

Arot[a_,b_,g_]:=rz[a].ry[b].rz[g]
Aroti[a_,b_,g_]:=rz[-g].ry[-b].rz[-a]

Define rotated shielding tensor

SA[sxx_,syy_,szz_,a_,b_,g_]:=Arot[a,b,g].(({{sxx, 0, 0}, {0, syy, 0}, {0, 0, szz}})).Aroti[a,b,g]

Define deviation of two matrices as the sum of the square of the differences between elements

chisqmat[mat1_,mat2_]:=Sum[j,{j,Sum[k,{k,(((mat1-mat2))^2)}]}]

Initial Steps:

Beginning with a Gaussian NMR calculation log file
	1.  Open the log file in GaussView to display the molecule 
	2.  Save the molecular coordinates in pdb format ("file", "save", "files of type", "pdb").
  	3.  From the Gaussian log file extract the Cartesian coordinates and the full chemical shielding tensor 
		for the atom(s) of interest.

Determine PAS components and Euler angles for orienting the tensor in the molecular/Gaussian frame

Notes :  Due to the symmetry of the shielding tensor, the Euler angles are not distinct.  Here the convention is chosen to find the solution between 0 ≤ α ≤ 2 π, 0 ≤ β ≤ π/2, and 0 ≤ γ ≤ π.

Modified to take the displacement tensor.

Isotropic=385;
XX=242; YX=0.0; ZX=0.0;
XY=0.0; YY=610; ZY=0.0;
XZ=0.0; YZ=0.0; ZZ=303;
atomcoord={0.000000, 0.000000, 0.000000};

ShieldingTensor=(({{XX, XY, XZ}, {YX, YY, YZ}, {ZX, ZY, ZZ}}));
ShieldingTensorSym=1/2(ShieldingTensor+Transpose[ShieldingTensor]);
({sxx,szz,syy}=Eigenvalues[ShieldingTensorSym])

{chisq,Euler}=NMinimize[{chisqmat[SA[sxx,syy,szz,a,b,g],ShieldingTensorSym],0≤a≤(2π)&&0≤b≤π/2&&0≤g≤π},{a,b,g},Method->"RandomSearch"]

Plot tensor in the molecular/Gaussian frame

Notes :  	1.  Set Shift = False 
		2.  const should be 0 for displacement tensors
		3.  scaling is used to adjust the size of the tensor relative to the molecule; 1 for displacement tensor
		4.  sa can be used to adjust the length of the principal axes; not used for displacement tensor

dir="C:\\temp\\";
structurefile="molecule.pdb";
pMolecule=Import[dir<>structurefile,"Rendering" -> "BallAndStick"]
prim=InputForm[%];
pMoleculewf=Import[dir<>structurefile,"Rendering" -> "Wireframe"]

scale=0.66;
pMoleculethin=prim[[1]]/.{Sphere[a_,b_]:>Sphere[a,scale b],Cylinder[a_,b_]:>Cylinder[a,scale b]}

Shift=False;
ref=0
const=0;
scaling=1;
sa=1.1;

{α,β,γ}={a,b,g}/.Euler;
{xc,yc,zc}=atomcoord 100;
If[Shift,{Sxx,Syy,Szz}=(ref-{sxx,syy,szz} + const)scaling,{Sxx,Syy,Szz}=({sxx,syy,szz} + const)scaling];

Modified to take the displacement tensor.

pShieldingTensor=Show[ContourPlot3D[((Aroti[α,β,γ].{x-xc,y-yc,z-zc})[[1]])^2/Sxx+((Aroti[α,β,γ].{x-xc,y-yc,z-zc})[[2]])^2/Syy+((Aroti[α,β,γ].{x-xc,y-yc,z-zc})[[3]])^2/Szz==1,{x,-300,300},{y,-300,300},{z,-300,300},Mesh->False,ContourStyle->{{Orange,Opacity[0.7],Specularity[White,20]}},PlotPoints->500],Background->White,Boxed->False,Axes->None,AxesStyle->White,Lighting->{{"Ambient",GrayLevel[.3]},{"Directional",White,ImageScaled[{1,1,1}]}},ImageSize->{500,357},PlotRange->{{-100,100},{-100,100},{-100,100}},SphericalRegion->True,(*ViewAngle->π/10,*)ImagePadding->{{0,0},{0,40}}(*,ViewPoint->{0,-7,-10}/10*)]

Show[pMoleculewf,pShieldingTensor,Background->White,Boxed->False,Axes->None,AxesStyle->White,Lighting->{{"Ambient",GrayLevel[.3]},{"Directional",White,ImageScaled[{1,1,1}]}},SphericalRegion->False,ImageSize->2{357,500},ViewPoint->{-1.6 ,2.6,-2.1}]

Show[pMoleculethin,pShieldingTensor,Background->White,Boxed->False,Axes->None,AxesStyle->White,Lighting->{{"Ambient",GrayLevel[.3]},{"Directional",White,ImageScaled[{1,1,1}]}},SphericalRegion->False,ImageSize->2{357,500},ViewPoint->{-1.6 ,2.6,-2.1}]