-
Notifications
You must be signed in to change notification settings - Fork 25
/
Copy pathexample_fromPDFgui.py
65 lines (52 loc) · 1.79 KB
/
example_fromPDFgui.py
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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
'''
This example will show how to perform mPDF refinements after having done the atomic
PDF refinement in PDFgui. We use data from MnO at 15 K.
'''
# Import necessary functions
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
from diffpy.mpdf import *
from diffpy.Structure import loadStructure
# Create the structure from our cif file, update the lattice params
structureFile = "MnO_R-3m.cif"
mnostructure = loadStructure(structureFile)
lat = mnostructure.lattice
lat.a,lat.b,lat.c = 3.1505626,3.1505626,7.5936979 ## refined values from PDFgui
# Create the Mn2+ magnetic species
mn2p = MagSpecies(struc=mnostructure, label='Mn2+', magIdxs=[0,1,2],
basisvecs=2.5*np.array([1,0,0]), kvecs=np.array([0,0,1.5]),
ffparamkey='Mn2')
# Create and prep the magnetic structure
mstr = MagStructure()
mstr.loadSpecies(mn2p)
mstr.makeAtoms()
mstr.makeSpins()
mstr.makeFF()
# Set up the mPDF calculator
mc = MPDFcalculator(magstruc=mstr, gaussPeakWidth=0.2)
# Load the data
PDFfitFile = 'MnOfit_PDFgui.fgr'
rexp,Drexp = getDiffData([PDFfitFile]) # this reads in the fit file
mc.rmin = rexp.min()
mc.rmax = rexp.max()
# Do the refinement
def residual(p, yexp, mcalc):
mcalc.paraScale, mcalc.ordScale = p
return yexp-mcalc.calc(both=True)[2]
p0 = [5.0,3.0] # initial parameter values (paraScale, ordScale)
pOpt = leastsq(residual, p0, args=(Drexp,mc))
print(pOpt)
fit=mc.calc(both=True)[2]
# Plot the results
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(rexp, Drexp, marker='o', mfc='none', mec='b', linestyle='none')
ax.plot(rexp, fit, 'r-', lw=2)
ax.set_xlim(xmin=mc.rmin, xmax=mc.rmax)
ax.set_xlabel(u'r (Å)')
ax.set_ylabel(u'd(r) (Å$^{-2}$)')
plt.show()