-
Notifications
You must be signed in to change notification settings - Fork 0
/
newHydronIsotopologue.m
157 lines (125 loc) · 5.55 KB
/
newHydronIsotopologue.m
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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
function Isotopologue = newHydronIsotopologue(Nuclei,System)
NuclearList = 1:Nuclei.number;
isSolvent = Nuclei.isSolvent(NuclearList);
switch System.HydrogenExchange
case 'OH'
for iNuc = NuclearList
if Nuclei.Exchangeable(iNuc)
deuteriumFraction = System.deuteriumFraction;
else
deuteriumFraction = System.deuteriumFraction_nonExchangeable;
end
type = Nuclei.Type{iNuc};
if isSolvent(iNuc) && (strcmp(type,'1H') || strcmp(type,'2H'))
if rand() > deuteriumFraction
Nuclei.Type{iNuc} = '1H';
Nuclei.Spin(iNuc) = 0.5; % hbar
Nuclei.StateMultiplicity(iNuc) = 2*Nuclei.Spin(iNuc) +1;
Nuclei.Nuclear_g(iNuc) = 5.58569;
Nuclei.NumberStates(iNuc) = int8(2);
% D =============================================================
else
Nuclei.Type{iNuc} = '2H';
Nuclei.Spin(iNuc) = 1; % hbar
Nuclei.StateMultiplicity(iNuc) = 2*Nuclei.Spin(iNuc) +1;
Nuclei.Nuclear_g(iNuc) = 0.857438;
Nuclei.NumberStates(iNuc) = int8(3);
end
if ~Nuclei.Exchangeable(iNuc)
sameMolecule = NuclearList(Nuclei.MoleculeID==Nuclei.MoleculeID(iNuc));
for jNuc = sameMolecule
if Nuclei.Exchangeable(jNuc) || iNuc==jNuc
continue;
end
Nuclei.Type{jNuc} = Nuclei.Type{iNuc};
Nuclei.Spin(jNuc) = Nuclei.Spin(iNuc); % hbar
Nuclei.StateMultiplicity(jNuc) = Nuclei.StateMultiplicity(iNuc);
Nuclei.Nuclear_g(jNuc) = Nuclei.Nuclear_g(iNuc);
Nuclei.NumberStates(jNuc) = Nuclei.NumberStates(iNuc);
end
end
end
end
case 'full'
for iNuc = 1:Nuclei.number
type = Nuclei.Type{iNuc};
if isSolvent(iNuc) && (strcmp(type,'1H') || strcmp(type,'2H'))
if rand() > System.deuteriumFraction
Nuclei.Type{iNuc} = '1H';
Nuclei.Spin(iNuc) = 0.5; % hbar
Nuclei.StateMultiplicity(iNuc) = 2*Nuclei.Spin(iNuc) +1;
Nuclei.Nuclear_g(iNuc) = 5.58569;
Nuclei.NumberStates(iNuc) = int8(2);
% D =============================================================
else
Nuclei.Type{iNuc} = '2H';
Nuclei.Spin(iNuc) = 1; % hbar
Nuclei.StateMultiplicity(iNuc) = 2*Nuclei.Spin(iNuc) +1;
Nuclei.Nuclear_g(iNuc) = 0.857438;
Nuclei.NumberStates(iNuc) = int8(3);
end
end
end
case 'none'
for iMol = Nuclei.MoleculeIDunique
MolecularNuclei = find(Nuclei.MoleculeID== iMol);
if isempty(MolecularNuclei)
continue;
end
testNuc = MolecularNuclei(1);
if isSolvent(testNuc)
for iNuc = MolecularNuclei
type = Nuclei.Type{iNuc};
if isSolvent(iNuc) && (strcmp(type,'1H') || strcmp(type,'2H'))
if rand() > System.deuteriumFraction
Nuclei.Type{iNuc} = '1H';
Nuclei.Spin(iNuc) = 0.5; % hbar
Nuclei.StateMultiplicity(iNuc) = 2*Nuclei.Spin(iNuc) +1;
Nuclei.Nuclear_g(iNuc) = 5.58569;
Nuclei.NumberStates(iNuc) = int8(2);
% D =============================================================
else
Nuclei.Type{iNuc} = '2H';
Nuclei.Spin(iNuc) = 1; % hbar
Nuclei.StateMultiplicity(iNuc) = 2*Nuclei.Spin(iNuc) +1;
Nuclei.Nuclear_g(iNuc) = 0.857438;
Nuclei.NumberStates(iNuc) = int8(3);
end
end
end
end
end
otherwise
error('Could not generate isotopologue.');
end
is1H = cellfun(@(x)isequal(x,'1H'),Nuclei.Type) & isSolvent;
is2H = cellfun(@(x)isequal(x,'2H'),Nuclei.Type) & isSolvent;
isExch = Nuclei.Exchangeable;
n1H = sum(is1H); % number of solvent protons.
n2H = sum(is2H); % number of solvent deuterons.
n1Hx = sum(is1H & isExch); % number of exchangeable solvent protons.
n2Hx = sum(is2H & isExch); % number of exchangeable solvent deuterons.
n1Hn = sum(is1H & ~isExch); % number of non-exchangeable solvent protons.
n2Hn = sum(is2H & ~isExch); % number of non-exchangeable solvent deuterons.
nH = (n1H+n2H); % number of solvent hydrons.
nHx = n1Hx + n2Hx; % number of exchangeable solvent hydrons.
nHn = n1Hn + n2Hn; % number of non-exchangeable solvent hydrons.
Ptarget = nHx/nH*System.deuteriumFraction + nHn/nH*System.deuteriumFraction_nonExchangeable;
P = n2H/(n1H+n2H);
Px = n2Hx/nHx;
Pn = n2Hn/nHn;
fprintf('Target: \n P(D) = %d/%d = %d. \n P(D|OH) = %d/%d = %d.\n P(D|CH) = %d/%d = %d.\n',...
round(nHx*System.deuteriumFraction + nHn*System.deuteriumFraction_nonExchangeable), nH, Ptarget, ...
round(nHx*System.deuteriumFraction), nHx, System.deuteriumFraction,...
round(nHn*System.deuteriumFraction_nonExchangeable), nHn, System.deuteriumFraction_nonExchangeable);
fprintf('Initialized %d hydrons: \n P(D) = %d/%d = %d. \n P(D|OH) = %d/%d = %d.\n P(D|CH) = %d/%d = %d.\n',nH,...
n2H, (n1H+n2H), P, ...
n2Hx, nHx, Px, ...
n2Hn, nHn, Pn);
Isotopologue = Nuclei;
Isotopologue.Isotopologue.Name = {'All','OH','CH'};
Isotopologue.Isotopologue.TypeNumber = [nH,nHx,nHn];
Isotopologue.Isotopologue.TargetFraction = [Ptarget,System.deuteriumFraction,System.deuteriumFraction_nonExchangeable];
Isotopologue.Isotopologue.Instance_2H_Number = [n2H,n2Hx,n2Hn];
Isotopologue.Isotopologue.InstanceFraction = [P,Px,Pn];
end