samerafach/Spintrum
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

First push, taken from the previous Spintrum software model

Browse Source
This commit is contained in:
Samer Afach
2017-01-10 11:14:51 +01:00
commit 7852c4d1d9
36 changed files with 168893 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

61
examples/ZULF/acetonitrile.py Normal file
View File

@@ -0,0 +1,61 @@
#!/usr/bin/python3
import spintrum
import math
import matplotlib.pyplot as plt
import numpy as np
import time
gammas = [4257.7e4,4257.7e4,4257.7e4,1070.8e4,1070.8e4,-431.6e4]
multips = [2,2,2,2,2,2]
jCouplings = \
[
[0,0,0,136.25,-9.94,-2.03],
[0,0,0,136.25,-9.94,-2.03],
[0,0,0,136.25,-9.94,-2.03],
[0,0,0,0,56.94,2.9],
[0,0,0,0,0,-17.53],
[0,0,0,0,0,0],
]
BThermal = 1.8e0
T2 = 10
sampleRate = 2e3
T= 1/sampleRate
points = 80000
gammah = 2*math.pi*4257.7
spinOp = spintrum.SpinOperations()
spinOp.add_operation(spintrum.SpinOperations.OPERATION__THERMAL_POPULATE,
{'Bx': 0, 'By': 0, 'Bz': BThermal, 'T': 293.778})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__TIP_SPINS,
{'direction': 'y', 'BVsTArea': 4*math.pi/gammah})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__SET_HAMILTONIAN,
{'Bx': 0, 'By': 0, 'Bz': 0})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION,
{'samplingRate': sampleRate, 'measurementDirection': 'z'})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT,
{'points': points, 'threads': 4})
start_time = time.time()
signal = spintrum.simulate(gyromagneticRatios=gammas,
jCouplings=jCouplings,
spinMultiplicities=multips,
spinOperations=spinOp)
signal = signal - np.mean(signal)
signal = [signal[i]*math.exp(-i/sampleRate/T2) for i in range(len(signal))]
print("Simulation lastet: " + repr(time.time()-start_time) + "s")
plt.plot(signal)
plt.show()
fft = spintrum.FFTSpectralDensity(signal, sampleRate)
plt.plot(fft['x'], fft['y'])
plt.show()

62
examples/ZULF/benzene.py Normal file
View File

@@ -0,0 +1,62 @@
#!/usr/bin/python3
import spintrum
import math
import matplotlib.pyplot as plt
import numpy as np
import time
gammas = [4257.7e4,4257.7e4,4257.7e4,4257.7e4,4257.7e4,4257.7e4,1070.8e4]
multips = [2,2,2,2,2,2,2]
jCouplings = \
[
[0,7.54,1.38,0.661,1.38,7.54,158.354],
[0,0,7.543,1.377,0.658,1.373,1.133],
[0,0,0,7.535,1.382,0.658,7.607],
[0,0,0,0,7.535,1.377,-1.296],
[0,0,0,0,0,7.543,7.607],
[0,0,0,0,0,0,1.133],
[0,0,0,0,0,0,0]
]
BThermal = 1.8e0
T2 = 10
sampleRate = 1e3
T= 1/sampleRate
points = 80000
gammah = 2*math.pi*4257.7
spinOp = spintrum.SpinOperations()
spinOp.add_operation(spintrum.SpinOperations.OPERATION__THERMAL_POPULATE,
{'Bx': 0, 'By': 0, 'Bz': BThermal, 'T': 293.778})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__TIP_SPINS,
{'direction': 'y', 'BVsTArea': 4*math.pi/gammah})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__SET_HAMILTONIAN,
{'Bx': 0, 'By': 0, 'Bz': 0})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION,
{'samplingRate': sampleRate, 'measurementDirection': 'z'})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT,
{'points': points, 'threads': 4})
start_time = time.time()
signal = spintrum.simulate(gyromagneticRatios=gammas,
jCouplings=jCouplings,
spinMultiplicities=multips,
spinOperations=spinOp)
signal = signal - np.mean(signal)
signal = [signal[i]*math.exp(-i/sampleRate/T2) for i in range(len(signal))]
print("Simulation lastet: " + repr(time.time()-start_time) + "s")
plt.plot(signal)
plt.show()
fft = spintrum.FFTSpectralDensity(signal, sampleRate)
plt.plot(fft['x'], fft['y'])
plt.show()

155
examples/ZULF/benzene_fit.py Normal file
View File

@@ -0,0 +1,155 @@
#!/usr/bin/python3
import spintrum
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
import multiprocessing
import mpmath
mpmath.mp.dps = 25
def filter_spectrum(spec, freq_lim):
x_axis = np.array([])
y_axis = np.array([])
for i in range(len(freq_lim)):
for j in range(len(spec["x"])):
if freq_lim[i][0] <= spec["x"][j] <= freq_lim[i][1]:
x_axis = np.append(x_axis,spec["x"][j])
y_axis = np.append(y_axis,spec["y"][j])
return {"x": x_axis, "y": y_axis}
with open("../data/benzeneSignal.txt") as f:
data = f.readlines()
data = np.array(list(map(np.double,data)))
data = data - np.mean(data)
#defining initial parameters
gammas = [4257.7e4,4257.7e4,4257.7e4,4257.7e4,4257.7e4,4257.7e4,1070.8e4]
multiplicities = [2, 2, 2, 2, 2, 2, 2]
gammah = 2*math.pi*4257.7e4
BThermal = 1.8e0
T2 = 1
points = len(data)
sample_rate = 2000
spectrum_range = [[5, 50], [80, 119], [121, 179], [181, 239], [241, 299]]
jCouplings = \
[
[0,0,0,0,0,0,0],
[0,0,0,0,0,0,0],
[0,0,0,0,0,0,0],
[0,0,0,0,0,0,0],
[0,0,0,0,0,0,0],
[0,0,0,0,0,0,0],
[0,0,0,0,0,0,0],
]
experimental_spectrum = filter_spectrum(spintrum.FFTSpectralDensity(data, sample_rate), spectrum_range)
gen = spintrum.SpinSimulator(gyromagneticRatios=gammas,
jCouplings=jCouplings,
spinMultiplicities=multiplicities,
doPrint=False)
spinOp = spintrum.SpinOperations()
spinOp.add_operation(spintrum.SpinOperations.OPERATION__THERMAL_POPULATE,
{'Bx': 0, 'By': 0, 'Bz': BThermal, 'T': 293.778})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__TIP_SPINS,
{'direction': 'y', 'BVsTArea': 4*math.pi/gammah})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__SET_HAMILTONIAN,
{'Bx': 0, 'By': 0, 'Bz': 0})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION,
{'samplingRate': sample_rate, 'measurementDirection': 'z'})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT,
{'points': points, 'threads': multiprocessing.cpu_count()})
gen.update_parameters(spinOperations=spinOp, jCouplings=jCouplings, gyromagneticRatios=gammas)
def generate_spectrum(params):
amplitude = params[13]
T2 = params[14]
jCouplings = \
[
[0, params[0], params[1], params[2],params[1],params[0],params[3]],
[0, 0, params[4], params[5], params[6], params[7], params[8]],
[0, 0, 0, params[9], params[10], params[6], params[11]],
[0, 0, 0, 0, params[9], params[5], params[12]],
[0, 0, 0, 0, 0, params[4], params[11]],
[0, 0, 0, 0, 0, 0, params[8]],
[0, 0, 0, 0, 0, 0, 0],
]
gen.update_parameters(jCouplings=jCouplings)
signal = amplitude*gen.simulate()
signal = signal - np.mean(signal)
signal = [signal[i] * math.exp(-i / sample_rate / T2) for i in range(len(signal))]
fft = filter_spectrum(spintrum.FFTSpectralDensity(signal, sample_rate), spectrum_range)
# plt.plot(fft['x'],fft['y'],realSpectrum['x'],realSpectrum['y'])
# plt.show()
return fft
def objective_func(params):
print(params)
spect = generate_spectrum(params)
func_value = np.sqrt(np.sum((spect['y'] - experimental_spectrum['y']) ** 2))
print('Objective function value:', func_value)
return func_value
pars = np.array([7.54,1.38,0.661,158.354,7.543,1.377,0.658,1.373,1.133,7.535,1.382,7.607,-1.296,0.00002,6])
#pars = np.array([-1.06785235e+00,2.26791767e+02,4.26684847e+00,4.53620118e-05,2.73284204e+00])
# get_objective_func(pars)
optimum_params = scipy.optimize.fmin(objective_func, pars, ftol = 5.0e-13, maxfun = 1)
opt_j_couplings = mpmath.matrix([
[0, optimum_params[0], optimum_params[1], optimum_params[2], optimum_params[1], optimum_params[0], optimum_params[3]],
[0, 0, optimum_params[4], optimum_params[5], optimum_params[6], optimum_params[7], optimum_params[8]],
[0, 0, 0, optimum_params[9], optimum_params[10], optimum_params[6], optimum_params[11]],
[0, 0, 0, 0, optimum_params[9], optimum_params[5], optimum_params[12]],
[0, 0, 0, 0, 0, optimum_params[4], optimum_params[11]],
[0, 0, 0, 0, 0, 0, optimum_params[8]],
[0, 0, 0, 0, 0, 0, 0],
])
CorrelationMatrix = spintrum.get_correlation_matrix(optimum_params, generate_spectrum, objective_func)
print('Correlation matrix: ', '\n', CorrelationMatrix)
print('Parameters at minimum: ','\n', optimum_params)
errorbars = spintrum.get_errorbars(CorrelationMatrix)
print('Standard deviations of parameters: ','\n',errorbars)
#writout to file
np.set_printoptions(suppress=True, precision=20)
fr = open("../data/benzeneFitResult.txt","w")
fr.write('Minimum at J-couplings:\n')
fr.write(str(opt_j_couplings) + '\n')
fr.write('Minimum at paramters:\n')
fr.write(str(optimum_params) + '\n')
fr.write('Parameters standard deviations:\n')
fr.write(str(errorbars)+'\n')
fr.write('Parameters correlation matrix:\n')
fr.write(str(CorrelationMatrix)+'\n')
fr.write('Minimum objective function value: ' + str(objective_func(optimum_params)))
fr.close()
fft_final = generate_spectrum(optimum_params)
plt.plot(fft_final['x'], fft_final['y'], experimental_spectrum['x'], experimental_spectrum['y'])
plt.show()

155
examples/ZULF/ethanol_fit.py Normal file
View File

@@ -0,0 +1,155 @@
#!/usr/bin/python3
import spintrum
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
import multiprocessing
import mpmath as mp
mp.mp.dps = 25
def filter_spectrum(spec, freq_lim):
x_axis = np.array([])
y_axis = np.array([])
for i in range(len(freq_lim)):
for j in range(len(spec["x"])):
if freq_lim[i][0] <= spec["x"][j] <= freq_lim[i][1]:
x_axis = np.append(x_axis,spec["x"][j])
y_axis = np.append(y_axis,spec["y"][j])
return {"x": x_axis, "y": y_axis}
#defining initial parameters
with open("../data/ethanolSignal.txt") as f:
data = f.readlines()
data = np.array(list(map(np.double,data)))
data = data - np.mean(data)
gammas = [4257.7e4,4257.7e4,4257.7e4,4257.7e4,4257.7e4,1070.8e4]
multips = [2,2,2,2,2,2]
gammah = 2*math.pi*4257.7e4
BThermal = 1.8e0
T2 = 1
points = len(data)
sampleRate = 2000
spectrumRange = [[0,299]]
jCouplings = \
[
[0,0,0,0,0,0],
[0,0,0,0,0,0],
[0,0,0,0,0,0],
[0,0,0,0,0,0],
[0,0,0,0,0,0],
[0,0,0,0,0,0],
]
experimental_spectrum = filter_spectrum(spintrum.FFTSpectralDensity(data, 2000), spectrumRange)
gen = spintrum.SpinSimulator(gyromagneticRatios=gammas,
jCouplings=jCouplings,
spinMultiplicities=multips,
doPrint=False)
spinOp = spintrum.SpinOperations()
spinOp.add_operation(spintrum.SpinOperations.OPERATION__THERMAL_POPULATE,
{'Bx': 0, 'By': 0, 'Bz': BThermal, 'T': 293.778})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__TIP_SPINS,
{'direction': 'y', 'BVsTArea': 4*math.pi/gammah})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__SET_HAMILTONIAN,
{'Bx': 0, 'By': 0, 'Bz': 0})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION,
{'samplingRate': sampleRate, 'measurementDirection': 'z'})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT,
{'points': points, 'threads': multiprocessing.cpu_count()})
gen.update_parameters(spinOperations=spinOp, jCouplings=jCouplings, gyromagneticRatios=gammas)
def generate_spectrum(params):
amplitude = params[3]
T2 = params[4]
jCouplings = \
[
[0, 0, 0, params[0],params[1],params[1]],
[0, 0, 0, params[0],params[1],params[1]],
[0, 0, 0, params[0],params[1],params[1]],
[0, 0, 0, 0, 0, params[2]],
[0, 0, 0, 0, 0, params[2]],
[0, 0, 0, 0, 0, 0],
]
gen.update_parameters(jCouplings=jCouplings)
signal = amplitude*gen.simulate()
signal = signal - np.mean(signal)
signal = [signal[i]*math.exp(-i/sampleRate/T2) for i in range(len(signal))]
fft = filter_spectrum(spintrum.FFTSpectralDensity(signal, 2000), spectrumRange)
# plt.plot(fft['x'],fft['y'],realSpectrum['x'],realSpectrum['y'])
# plt.show()
return fft
def objective_func(params):
print(params)
spect = generate_spectrum(params)
funcValue = np.sqrt(np.sum((spect['y'] - experimental_spectrum['y']) ** 2))
print('Objective function value: ',funcValue)
return funcValue
pars = np.array([7.1,141.0,-4.65,0.00002,6])
#pars = np.array([-1.06785235e+00,2.26791767e+02,4.26684847e+00,l4.53620118e-05,2.73284204e+00])
# get_objective_func(pars)
optimum_params = scipy.optimize.fmin(objective_func, pars, ftol = 5.0e-13, maxfun = 1)
opt_j_couplings = mp.matrix([
[0, 0, 0, optimum_params[0], optimum_params[1], optimum_params[1]],
[0, 0, 0, optimum_params[0], optimum_params[1], optimum_params[1]],
[0, 0, 0, optimum_params[0], optimum_params[1], optimum_params[1]],
[0, 0, 0, 0, 0, optimum_params[2]],
[0, 0, 0, 0, 0, optimum_params[2]],
[0, 0, 0, 0, 0, 0],
])
CorrelationMatrix = spintrum.get_correlation_matrix(optimum_params, generate_spectrum, objective_func)
print('Correlation matrix: ','\n',CorrelationMatrix)
print('Parameters at minimum: ','\n', optimum_params)
errorbars = spintrum.get_errorbars(CorrelationMatrix)
print('Standard deviations of parameters: ','\n',errorbars)
#writout to file
np.set_printoptions(suppress=True,precision=20)
fr = open("../data/ethanolFitResult.txt","w")
fr.write('Minimum at J-couplings:\n')
fr.write(str(opt_j_couplings) + '\n')
fr.write('Minimum at paramters:\n')
fr.write(str(optimum_params) + '\n')
fr.write('Parameters standard deviations:\n')
fr.write(str(errorbars)+'\n')
fr.write('Parameters correlation matrix:\n')
fr.write(str(CorrelationMatrix)+'\n')
fr.write('Minimum objective function value: ' + str(objective_func(optimum_params)))
fr.close()
fft_final = generate_spectrum(optimum_params)
plt.plot(fft_final['x'], fft_final['y'], experimental_spectrum['x'], experimental_spectrum['y'])
plt.show()

80
examples/ZULF/formamide.py Normal file
View File

@@ -0,0 +1,80 @@
#!/usr/bin/python3
import spintrum
import math
import matplotlib.pyplot as plt
import numpy as np
import time
gammas = [-431.6e4,4257.7e4,4257.7e4,4257.7e4]
multips = [2,2,2,2]
jCouplings = \
[
[0,-89.3,-89.3,-13.5],
[0,0,0,8],
[0,0,0,8],
[0,0,0,0],
]
BThermal = 1.8e0
T2 = 5
sampleRate = 2e3
T= 1/sampleRate
points = 80000
specRange = [[0,50],[100,200]]
gammah = 2*math.pi*4257.7e4
spinOp = spintrum.SpinOperations()
spinOp.add_operation(spintrum.SpinOperations.OPERATION__THERMAL_POPULATE,
{'Bx': 0, 'By': 0, 'Bz': BThermal, 'T': 293.778})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__TIP_SPINS,
{'direction': 'y', 'BVsTArea': 4*math.pi/gammah})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__SET_HAMILTONIAN,
{'Bx': 0, 'By': 0, 'Bz': 0})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION,
{'samplingRate': sampleRate, 'measurementDirection': 'z'})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT,
{'points': points, 'threads': 4})
def filterSpectrum(spec, freqLim):
xAxis = np.array([])
yAxis = np.array([])
for i in range(len(freqLim)):
for j in range(len(spec["x"])):
if spec["x"][j] >= freqLim[i][0] and spec["x"][j] <= freqLim[i][1]:
xAxis = np.insert(xAxis,-0,spec["x"][j])
yAxis = np.insert(yAxis,-0,spec["y"][j])
return {"x": xAxis, "y": yAxis}
start_time = time.time()
signal = spintrum.simulate(gyromagneticRatios=gammas,
jCouplings=jCouplings,
spinMultiplicities=multips,
spinOperations=spinOp)
print("Simulation lastet: " + repr(time.time()-start_time) + "s")
signal = signal - np.mean(signal)
signal = [signal[i]*math.exp(-i/sampleRate/T2) for i in range(len(signal))]
plt.plot(signal)
plt.show()
FFT = spintrum.FFTSpectralDensity(signal, sampleRate)
fft = filterSpectrum(FFT,specRange)
plt.plot(fft['x'], fft['y'])
plt.show()

142
examples/ZULF/formicacid_fit.py Normal file
View File

@@ -0,0 +1,142 @@
#!/usr/bin/python3
import spintrum
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
import multiprocessing
import mpmath as mp
def filter_spectrum(spec, freq_lim):
x_axis = np.array([])
y_axis = np.array([])
for i in range(len(freq_lim)):
for j in range(len(spec["x"])):
if freq_lim[i][0] <= spec["x"][j] <= freq_lim[i][1]:
x_axis = np.append(x_axis,spec["x"][j])
y_axis = np.append(y_axis,spec["y"][j])
return {"x": x_axis, "y": y_axis}
with open("../data/formicacidSignal.txt") as f:
data = f.readlines()
data = np.array(list(map(np.double,data)))
data = data - np.mean(data)
gammas = [4257.7e4,1070.8e4]
multips = [2,2]
gammah = 2*math.pi*4257.7e4
BThermal = 1.8e0
T2 = 1
points = len(data)
sampleRate = 2000
spectrumRange = [[210,230]]
jCouplings = \
[
[0,0],
[0,0],
]
realSpectrum = filter_spectrum(spintrum.FFTSpectralDensity(data, 2000), spectrumRange)
gen = spintrum.SpinSimulator(gyromagneticRatios=gammas,
jCouplings=jCouplings,
spinMultiplicities=multips,
doPrint=False)
spinOp = spintrum.SpinOperations()
spinOp.add_operation(spintrum.SpinOperations.OPERATION__THERMAL_POPULATE,
{'Bx': 0, 'By': 0, 'Bz': BThermal, 'T': 293.778})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__TIP_SPINS,
{'direction': 'y', 'BVsTArea': 4*math.pi/gammah})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__SET_HAMILTONIAN,
{'Bx': 0, 'By': 0, 'Bz': 0})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION,
{'samplingRate': sampleRate, 'measurementDirection': 'z'})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT,
{'points': points, 'threads': multiprocessing.cpu_count()})
gen.update_parameters(spinOperations=spinOp, jCouplings=jCouplings, gyromagneticRatios=gammas)
def generate_spectrum(params):
amplitude = params[1]
T2 = params[2]
jCouplings = \
[
[0, params[0]],
[0, 0],
]
gen.update_parameters(jCouplings=jCouplings)
signal = amplitude*gen.simulate()
signal = signal - np.mean(signal)
signal = [signal[i]*math.exp(-i/sampleRate/T2) for i in range(len(signal))]
fft = filter_spectrum(spintrum.FFTSpectralDensity(signal, 2000), spectrumRange)
# plt.plot(fft['x'],fft['y'],realSpectrum['x'],realSpectrum['y'])
# plt.show()
return fft
def objective_func(params):
print(params)
spect = generate_spectrum(params)
funcValue = np.sqrt(np.sum((spect['y'] - realSpectrum['y'])**2))
print('Objective function value: ',funcValue)
return funcValue
pars = np.array([222.5,0.000003,6])
#pars = np.array([-1.06785235e+00,2.26791767e+02,4.26684847e+00,4.53620118e-05,2.73284204e+00])
# get_objective_func(pars)
optimum_params = scipy.optimize.fmin(objective_func, pars, ftol = 5.0e-13, maxfun = 1e6)
optjCouplings = mp.matrix([
[0, optimum_params[0]],
[0, 0],
])
CorrelationMatrix = spintrum.get_correlation_matrix(optimum_params, generate_spectrum, objective_func)
print('Correlation matrix: ','\n',CorrelationMatrix)
print('Parameters at minimum: ','\n', optimum_params)
errorbars = spintrum.get_errorbars(CorrelationMatrix)
print('Standard deviations of parameters: ','\n',errorbars)
#writout to file
np.set_printoptions(suppress=True,precision=20)
fr = open("../data/formicacidFitResult.txt","w")
fr.write('Minimum at J-couplings:\n')
fr.write(str(optjCouplings)+'\n')
fr.write('Minimum at paramters:\n')
fr.write(str(optimum_params) + '\n')
fr.write('Parameters standard deviations:\n')
fr.write(str(errorbars)+'\n')
fr.write('Parameters correlation matrix:\n')
fr.write(str(CorrelationMatrix)+'\n')
fr.write('Minimum objective function value: ' + str(objective_func(optimum_params)))
fr.close()
fft_final = generate_spectrum(optimum_params)
plt.plot(fft_final['x'],fft_final['y'],realSpectrum['x'],realSpectrum['y'])
plt.show()

60
examples/ZULF/methylformate.py Normal file
View File

@@ -0,0 +1,60 @@
#!/usr/bin/python3
import spintrum
import math
import matplotlib.pyplot as plt
import numpy as np
import time
gammas = [4257.7e4,4257.7e4,4257.7e4,4257.7e4,1070.8e4]
multips = [2,2,2,2,2]
jCouplings = \
[
[0,-0.8,-0.8,-0.8,226.8],
[0,0,0,0,4.0],
[0,0,0,0,4.0],
[0,0,0,0,4.0],
[0,0,0,0,0],
]
BThermal = 1.8e0
T2 = 10
sampleRate = 2e3
T= 1/sampleRate
points = 80000
gammah = 2*math.pi*4257.7
spinOp = spintrum.SpinOperations()
spinOp.add_operation(spintrum.SpinOperations.OPERATION__THERMAL_POPULATE,
{'Bx': 0, 'By': 0, 'Bz': BThermal, 'T': 293.778})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__TIP_SPINS,
{'direction': 'y', 'BVsTArea': 4*math.pi/gammah})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__SET_HAMILTONIAN,
{'Bx': 0, 'By': 0, 'Bz': 0})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION,
{'samplingRate': sampleRate, 'measurementDirection': 'z'})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT,
{'points': points, 'threads': 4})
start_time = time.time()
signal = spintrum.simulate(gyromagneticRatios=gammas,
jCouplings=jCouplings,
spinMultiplicities=multips,
spinOperations=spinOp)
signal = signal - np.mean(signal)
signal = [signal[i]*math.exp(-i/sampleRate/T2) for i in range(len(signal))]
print("Simulation lastet: " + repr(time.time()-start_time) + "s")
plt.plot(signal)
plt.show()
fft = spintrum.FFTSpectralDensity(signal, sampleRate)
plt.plot(fft['x'], fft['y'])
plt.show()

156
examples/ZULF/methylformate_fit.py Normal file
View File

@@ -0,0 +1,156 @@
#!/usr/bin/python3
import spintrum
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
import multiprocessing
import mpmath as mp
def filter_spectrum(spec, freq_lim):
x_axis = np.array([])
y_axis = np.array([])
for i in range(len(freq_lim)):
for j in range(len(spec["x"])):
if freq_lim[i][0] <= spec["x"][j] <= freq_lim[i][1]:
x_axis = np.append(x_axis,spec["x"][j])
y_axis = np.append(y_axis,spec["y"][j])
return {"x": x_axis, "y": y_axis}
#defining initial parameters
with open("../data/methylformateSignal.txt") as f:
data = f.readlines()
data = np.array(list(map(np.double,data)))
data = data - np.mean(data)
gammas = [4257.7e4,4257.7e4,4257.7e4,4257.7e4,1070.8e4]
multips = [2,2,2,2,2]
gammah = 2*math.pi*4257.7e4
BThermal = 1.8e0
T2 = 1
points = len(data)
sampleRate = 2000
spectrumRange = [[70,119],[121,179],[181,239],[241,299]]
jCouplings = \
[
[0,0,0,0,0],
[0,0,0,0,0],
[0,0,0,0,0],
[0,0,0,0,0],
[0,0,0,0,0],
]
experimental_spectrum = filter_spectrum(spintrum.FFTSpectralDensity(data, 2000), spectrumRange)
gen = spintrum.SpinSimulator(gyromagneticRatios=gammas,
jCouplings=jCouplings,
spinMultiplicities=multips,
doPrint=False)
spinOp = spintrum.SpinOperations()
spinOp.add_operation(spintrum.SpinOperations.OPERATION__THERMAL_POPULATE,
{'Bx': 0, 'By': 0, 'Bz': BThermal, 'T': 293.778})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__TIP_SPINS,
{'direction': 'y', 'BVsTArea': 4*math.pi/gammah})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__SET_HAMILTONIAN,
{'Bx': 0, 'By': 0, 'Bz': 0})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION,
{'samplingRate': sampleRate, 'measurementDirection': 'z'})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT,
{'points': points, 'threads': multiprocessing.cpu_count()})
gen.update_parameters(spinOperations=spinOp, jCouplings=jCouplings, gyromagneticRatios=gammas)
def generate_spectrum(params):
amplitude = params[3]
T2 = params[4]
jCouplings = \
[
[0, params[0], params[0], params[0],params[1]],
[0, 0, 0, 0, params[2]],
[0, 0, 0, 0, params[2]],
[0, 0, 0, 0, params[2]],
[0, 0, 0, 0, 0],
]
spinOp = spintrum.SpinOperations()
spinOp.add_operation(spintrum.SpinOperations.OPERATION__THERMAL_POPULATE,
{'Bx': 0, 'By': 0, 'Bz': BThermal, 'T': 293.778})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__TIP_SPINS,
{'direction': 'y', 'BVsTArea': 4 * math.pi / gammah})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__SET_HAMILTONIAN,
{'Bx': np.float(params[5]), 'By': 0, 'Bz': np.float(params[6])})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION,
{'samplingRate': sampleRate, 'measurementDirection': 'z'})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT,
{'points': points, 'threads': multiprocessing.cpu_count()})
gen.update_parameters(jCouplings=jCouplings, spinOperations=spinOp)
signal = amplitude*gen.simulate()
signal = signal - np.mean(signal)
signal = [signal[i]*math.exp(-i/sampleRate/T2) for i in range(len(signal))]
fft = filter_spectrum(spintrum.FFTSpectralDensity(signal, 2000), spectrumRange)
# plt.plot(fft['x'],fft['y'],realSpectrum['x'],realSpectrum['y'])
# plt.show()
return fft
def objective_func(params):
print(params)
spect = generate_spectrum(params)
funcValue = np.sqrt(np.sum((spect['y'] - experimental_spectrum['y']) ** 2))
print('Objective function value: ',funcValue)
return funcValue
pars = np.array([-1.06785235e+00,2.26791767e+02,4.26684847e+00,4.53620118e-05,2.73284204e+00,3.11541220e-09,2.80718054e-09])
#pars = np.array([-1.06785235e+00,2.26791767e+02,4.26684847e+00,4.53620118e-05,2.73284204e+00])
# get_objective_func(pars)
optimum_params = scipy.optimize.fmin(objective_func, pars, ftol = 5.0e-13, maxfun = 1)
optjCouplings = mp.matrix([[0, optimum_params[0], optimum_params[0], optimum_params[0], optimum_params[1]], [0, 0, 0, 0, optimum_params[2]], [0, 0, 0, 0, optimum_params[2]], [0, 0, 0, 0, optimum_params[2]], [0, 0, 0, 0, 0]])
CorrelationMatrix = spintrum.get_correlation_matrix(optimum_params, generate_spectrum, objective_func)
print('Correlation matrix: ','\n',CorrelationMatrix)
print('Parameters at minimum: ','\n', optimum_params)
errorbars = spintrum.get_errorbars(CorrelationMatrix)
print('Standard deviations of parameters: ','\n',errorbars)
#writout to file
np.set_printoptions(suppress=True,precision=20)
fr = open("../data/methylformateFitResult.txt","w")
fr.write('Minimum at J-couplings:\n')
fr.write(str(optjCouplings)+'\n')
fr.write('Minimum at paramters:\n')
fr.write(str(optimum_params) + '\n')
fr.write('Parameters standard deviations:\n')
fr.write(str(errorbars)+'\n')
fr.write('Parameters correlation matrix:\n')
fr.write(str(CorrelationMatrix)+'\n')
fr.write('Minimum objective function value: ' + str(objective_func(optimum_params)))
fr.close()
fft_final = generate_spectrum(optimum_params)
plt.plot(fft_final['x'], fft_final['y'], experimental_spectrum['x'], experimental_spectrum['y'])
plt.show()

64
examples/ZULF/toluene.py Normal file
View File

@@ -0,0 +1,64 @@
#!/usr/bin/python3
import spintrum
import math
import matplotlib.pyplot as plt
import numpy as np
import time
gammas = [1070.8e4,4257.7e4,4257.7e4,4257.7e4,4257.7e4,4257.7e4,4257.7e4,4257.7e4,4257.7e4]
multips = [2,2,2,2,2,2,2,2,2]
jCouplings = \
[
[0,125.99,125.99,125.99,4.53,0.63,0.56,0.63,4.53],
[0,0,0,0,-0.69,0.3,-0.52,0.3,-0.69],
[0,0,0,0,-0.69,0.3,-0.52,0.3,-0.69],
[0,0,0,0,-0.69,0.3,-0.52,0.3,-0.69],
[0,0,0,0,0,7.655,1.273,0.61,1.902],
[0,0,0,0,0,0,7.417,1.442,0.61],
[0,0,0,0,0,0,0,7.417,1.273],
[0,0,0,0,0,0,0,0,7.655],
[0,0,0,0,0,0,0,0,0],
]
BThermal = 1.8e0
T2 = 10
sampleRate = 2e3
T= 1/sampleRate
points = 80000
gammah = 2*math.pi*4257.7
spinOp = spintrum.SpinOperations()
spinOp.add_operation(spintrum.SpinOperations.OPERATION__THERMAL_POPULATE,
{'Bx': 0, 'By': 0, 'Bz': BThermal, 'T': 293.778})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__TIP_SPINS,
{'direction': 'y', 'BVsTArea': 4*math.pi/gammah})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__SET_HAMILTONIAN,
{'Bx': 0, 'By': 0, 'Bz': 0})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION,
{'samplingRate': sampleRate, 'measurementDirection': 'z'})
spinOp.add_operation(spintrum.SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT,
{'points': points, 'threads': 4})
start_time = time.time()
signal = spintrum.simulate(gyromagneticRatios=gammas,
jCouplings=jCouplings,
spinMultiplicities=multips,
spinOperations=spinOp)
signal = signal - np.mean(signal)
signal = [signal[i]*math.exp(-i/sampleRate/T2) for i in range(len(signal))]
print("Simulation lastet: " + repr(time.time()-start_time) + "s")
#plt.plot(signal)
#plt.show()
fft = spintrum.FFTSpectralDensity(signal, sampleRate)
#plt.plot(fft['x'], fft['y'])
#plt.show()