CluE-MATLAB/setUpSystem.m at main · StollLab/CluE-MATLAB · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

function [System, Tensors, Nuclei,Clusters] = setUpSystem(System,Data)

tic
Method = [];
% set defaults base on specified parameters and for unspecified parameters
[System, Method, Data] = setDefaults(System,Method,Data);

if isfield(System,'Alpha')
 Alpha =System.Alpha;
else
 Alpha =0;
end
if isfield(System,'Beta')
 Beta =System.Gamma;
else
 Beta =0;
end
if isfield(System,'Gamma')
 Gamma =System.Gamma;
else
 Gamma =0;
end

% Get input data.
if ~isfield(Data,'InputData') || ~ischar(Data.InputData)
 error(['Data.InputData is not a valid file identifier.']);
end
InputData =Data.InputData;
if ~isfile(InputData)
 error(['Could not find the specified input file ', InputData, '.']);
end


% Set verbosity.
verbose =Method.verbose;

ifverbose, fprintf('Setting up...\n'); end

if min( (InputData(end-3:end)) =='.pdb')
 Nuclei = parseNuclei(System, Method, InputData);
 System.Electron.Coordinates = [0,0,0];
else
 error('Input data not recognized')
end
Method.order = min(Method.order,Nuclei.number);

ifverbose, fprintf('Setup initialized %i nuclei.\n', Nuclei.number); end

% Get rotation matrix from PDB frame to lab frame, via Euler angles
R_pdb2lab = rotateZYZ(Alpha,Beta,Gamma);

% Rotate nuclear coordinates.
Nuclei.Coordinates =Nuclei.Coordinates*R_pdb2lab';

% Rotate nuclear quadrupole tensors.
for inucleus =1:Nuclei.number
Nuclei.Qtensor(:,:,inucleus) =R_pdb2lab*Nuclei.Qtensor(:,:,inucleus)*R_pdb2lab';
% Elementwise Qtensor manipulation used for testing. The default filer is ones(3);
Nuclei.Qtensor(:,:,inucleus) =Nuclei.Qtensor(:,:,inucleus).*System.nuclear_quadrupole_filter;
end

% Rotate the g-matrix.
if isfield(System,'gFrame')
% Get rotation matrix.
 g2MolRotation = rotateZYZ(-System.gFrame(1),-System.gFrame(2),-System.gFrame(3));

% Rotate to molecular frame.
 System.gMatrix =g2MolRotation*System.gMatrix_gFrame*g2MolRotation';
end

% Rotate the g-matrix to the lab frame.
System.gMatrix =R_pdb2lab*System.gMatrix_gFrame*R_pdb2lab';

geff=System.gMatrix(3,3);

Tensors = pairwisetensors(Nuclei.Nuclear_g, Nuclei.Coordinates,...
 [1:Nuclei.number],Nuclei.Atensor, System.magneticField, System.ge, geff, System.muB, System.muN, System.mu0, System.hbar,System.theory,[]);
% ========================================================================
% Compile list of connected clusters
% ========================================================================
Clusters = [];
% Loop over cluster sizes, start at the largest (most time consuming) size
for clusterSize =Method.order:-1:1


ifverbose, fprintf('Finding clusters of size %d.\n', clusterSize); end

 Clusters{clusterSize} = findClusters(Nuclei, clusterSize);

% save the number of clusters of each size for export to the user
Nuclei.numberClusters(clusterSize) = size(Clusters{clusterSize},1);

ifverbose
 fprintf(' Found %d clusters of size %d.\n', size(Clusters{clusterSize},1),clusterSize);
end

end


% Determine the output file name.
if isfield(Data,'OutputData') && ~isempty(Data.OutputData)
 OutputData =Data.OutputData;
if ~strcmp(OutputData(end-3:end),'.mat')
 OutputData = [OutputData, '.mat'];
end
else
 OutputData ='';
end

% Save if filename is provided.
if ~isempty(OutputData)
 save(OutputData,'System','Tensors','Nuclei','Clusters','InputData');
end
toc
end