Semi-empirical mass formula fit

  • Difficulty: 028/100
  • Author: jeremy theler
  • Keywords: FIT, TO, VIA, VERBOSE, SEP, TEXT, FILE,

In 2010 I went to see a show by Sir Paul McCartney. One of the many things that moved me to tears was the fact that man in his seventies is still consistently doing and thinking about the same stuff he used to do and to think about fifty years ago. And this report that illustrates one of the applications of “wasora”: “http://www.talador.com.ar/jeremy/wasora” is the same stuff I have been thinking about more than ten years ago. This time, I got it right (I think).

During my second year in college I was interested in fitting functions to experimental data, in particular polynomials. I developed a method with the corresponding “code”: “http://www.talador.com.ar/jeremy/files/pfit-0.1.tar.gz” that now I can see was rather lame but at that time seemed a great discovery. When I switched to my definite “alma matter”: “http://www.ib.edu.ar”, I finally understood the proper way to fit a polynomial to scattered data. By that time, when studying nuclear physics, I also stumbled upon the “semi-empirical mass formula” :“http://en.wikipedia.org/wiki/Semi-empirical_mass_formula”. It was then time to combine these two brand-new concepts and try to obtain the damn coefficients writing another “lame academic paper”: “http://www.talador.com.ar/jeremy/files/yasmf.pdf”.

Ten years later, it is time to get it done the way it was supposed to from the very beginning: perform a non-linear multidimensional fit using “wasora”: “http://www.talador.com.ar/jeremy/wasora”.

fit.was

Given a certain nuclide characterized by the number~\(Z\) of protons and the total number~\(A\) of nucleons (i.e. protons plus neutrons), the “Weizsäcker’s mass formula” :“http://en.wikipedia.org/wiki/Semi-empirical_mass_formula” based on the liquid drop model gives a expected dependence of the binding energy per nucleon~\(B/A\) in terms of algebraic sums of powers of~\(Z\) and~\(A\) that involve some multiplicative factors that are to be determined experimentally.

FIGURE: [028-mass/fsm] A graphical representation of the binding energy per nucleon as a function of the mass number of the nuclide using the liquid drop model. This curve actually illustrates why both fission of heavy nuclei and fusion of light nuclei release energy.

To make a long story short, the semi-empirical mass formula says that the binding energy per nucleon of a nuclide of mass number~\(A\) with atomic number~\(Z\) is

!bt [ (A,Z) a_1 - a_2 A^{-1/3} - a_3 Z^2 A^{-4/3} - a_4 (A-2Z)^2 A^{-2} + a_5 A^{-} ] !et

with

!bt [ = \begin{cases} +1 & \text{for even-$A$ and even-$Z$} \\ 0 & \text{for odd-$A$} \\ -1 & \text{for even-$A$ and odd-$Z$} \\ \end{cases}

] !et

and~\(a_i\) for~\(i=1,\dots,5\) and~\(\gamma\) are six real constants to be empirically determined.

In this case, measured data from the “2003 Atomic Mass Evaluation”: “http://www.nndc.bnl.gov/masses/” is used as the experimental data. A point-wise defined function~\(D(A,Z)\) is defined, reading the data from a file with three columns, namely \(A\), \(Z\) and \(B/A\) called binding_per_A.dat. Another algebraic function~\(W(A,Z)\) is defined following the proposed functional form and leaving the to-be-determined constants as parameters. Using the FIT keyword, we ask “wasora”: “http://www.talador.com.ar/jeremy/wasora” (actually “GSL”: “http://www.gnu.org/software/gsl/”) to find the values of the constants that better fit \(W(A,Z)\) to \(D(A,Z)\). Using the VERBOSE keyword we can see how the iterative procedure advances. When a result is finally obtained, we write a file called the14.was that we will use later to compare our result with other fits.

FIGURE: [028-mass/experimental2d] Experimental binding energy for 2011 nuclides. Data taken from the “2003 Atomic Mass Evaluation”: “http://www.nndc.bnl.gov/masses/”

FIGURE: [028-mass/experimental3d] Surface view of experimental binding energy for 2011 nuclides. Data taken from the “2003 Atomic Mass Evaluation”: “http://www.nndc.bnl.gov/masses/”

# initial guess
a1 = 1
a2 = 1
a3 = 1
a4 = 1
a5 = 1
gamma = 1.5

# the functional form of weiszacker's formula
delta(A,Z) := if(is_odd(A), 0, if(is_even(Z), +1, -1))
W(A,Z) := a1 - a2*A^(-1/3) - a3*Z*(Z-1)*A^(-4/3) - a4*(A-2*Z)^2*A^(-2) + delta(A,Z) * a5*A^(-gamma)

# the experimental data
FUNCTION D(A,Z) FILE_PATH binding_per_A.dat

# fit W to D using the six parameters (and be verbose!)
FIT W TO D VIA a1 a2 a3 a4 a5 gamma VERBOSE

# once we are done, we write into the14.was what we found
# so we can afterward include this file to evaluate W(A,Z)
IF done_outer
 OUTPUT_FILE out the14.was
 PRINT SEP " " FILE out TEXT "a1 =" %g a1
 PRINT SEP " " FILE out TEXT "a2 =" %g a2
 PRINT SEP " " FILE out TEXT "a3 =" %g a3
 PRINT SEP " " FILE out TEXT "a4 =" %g a4
 PRINT SEP " " FILE out TEXT "a5 =" %g a5
 PRINT SEP " " FILE out TEXT "gamma =" %g  gamma
ENDIF
$ head binding_per_A_gr15.dat
head: cannot open ‘binding_per_A_gr15.dat’ for reading: No such file or directory
$ tail binding_per_A_gr15.dat
tail: cannot open ‘binding_per_A_gr15.dat’ for reading: No such file or directory
$ wasora fit.was
fititeration: 0  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.500e+00   status = success    |r|^2 = 526.8
fititeration: 1  1.407e+01  1.383e+01  5.787e-01  1.604e+01  4.051e+00  1.742e-01   status = success    |r|^2 = 57.1764
fititeration: 2  1.417e+01  1.409e+01  5.879e-01  1.612e+01  6.721e-01  2.076e-01   status = success    |r|^2 = 9.44812
fititeration: 3  1.417e+01  1.408e+01  5.875e-01  1.611e+01  7.540e-01  4.342e-01   status = success    |r|^2 = 6.36797
fititeration: 4  1.413e+01  1.399e+01  5.844e-01  1.608e+01  1.478e+00  8.361e-01   status = success    |r|^2 = 5.7591
fititeration: 5  1.409e+01  1.387e+01  5.801e-01  1.604e+01  3.153e+00  1.239e+00   status = success    |r|^2 = 5.67776
fititeration: 6  1.407e+01  1.383e+01  5.789e-01  1.605e+01  4.138e+00  1.261e+00   status = success    |r|^2 = 5.65002
fititeration: 7  1.407e+01  1.383e+01  5.789e-01  1.605e+01  4.152e+00  1.258e+00   status = success    |r|^2 = 5.64999
fititeration: 8  1.407e+01  1.383e+01  5.789e-01  1.605e+01  4.151e+00  1.258e+00   status = success    |r|^2 = 5.64999
fititeration: 9  1.407e+01  1.383e+01  5.789e-01  1.605e+01  4.151e+00  1.258e+00   status = success    |r|^2 = 5.64999
fititeration: 10     1.407e+01  1.383e+01  5.789e-01  1.605e+01  4.151e+00  1.258e+00   status = success    |r|^2 = 5.64999
# chisq/dof = 0.0159214
# a1    = 14.073 # +/- 0.27312
# a2    = 13.831 # +/- 0.69449
# a3    = 0.57892 # +/- 0.030422
# a4    = 16.046 # +/- 0.83499
# a5    = 4.1506 # +/- 3.2044
# gamma = 1.2578 # +/- 0.45427
$ 

Note that Weizsäcker’s formula is strongly non-linear in the spin term. Nevertheless, the fit managed to arrive at a solution that was relatively far away from the initial guess. The units of the constants \(a_i\) are~MeV (the same than the units of the third column of the file binding_per_A.dat). The reported uncertainties are to be taken as a relative measure of how pertinent are each of the fitted parameters for the experimental data set.

compare.was

We would like to compare first how well our fit resembles the experimental data, but also how other people’s fits do. In order to do this, we prepare one input file for each set of coefficients we want to compare. In particular, the file the14.was was automatically generated by the fit we just did:

@@@CODE 028-mass/the14.was envir=wasora

Emilio Segre in his book “Nuclei and particles” of 1965 reports (seg65.was):

@@@CODE 028-mass/seg65.was envir=wasora

Bernard L. Cohen in his book “Concepts of nuclear physic” of 1971 reports (coh71.was):

@@@CODE 028-mass/coh71.was envir=wasora

J. M. Pearson in his book “Nuclear physics: energy and matter” of 1986 reports (pea86.was):

@@@CODE 028-mass/pea86.was envir=wasora

In my undergraduate years using a lame fitting method, I used to report (the05.was):

@@@CODE 028-mass/the05.was envir=wasora

The following input called compare.was defines again the function \(W(A,Z)\) and uses one of the sets of coefficients listed above according to the parameter given in the command line. It also defines again the two-dimensional point-wise function \(D(A,Z)\) with the experimental data. Then it prints to an output file—whose name is taken from the command line—the experimental data, the estimation given by the semi-empirical mass formula with the chosen parameters and the absolute difference with respect to the experimental data. The first function that appears in the PRINT_FUNCTION instruction is \(D(A,Z)\) which is pointwise-defined so there is no need to give an explicit range. Therefore, both~\(W(A,Z)\) and~\(D(A,Z)-W(A,Z)\) are shown only at the points \((A,Z)\) for which experimental data exists.

# read the coefficients asked for in the command line
INCLUDE $1.was

# weiszacker's formula
delta(A,Z) := if(is_odd(A), 0, if(is_even(Z), +1, -1))
W(A,Z) := a1 - a2*A^(-1/3) - a3*Z*(Z-1)*A^(-4/3) - a4*(A-2*Z)^2*A^(-2) + delta(A,Z) * a5*A^(-gamma)

# experimental data
FUNCTION D(A,Z) FILE_PATH binding_per_A.dat

# write the functions to a file (D comes first so no explicit range is needed)
PRINT_FUNCTION D W D(A,Z)-W(A,Z) (D(A,Z)-W(A,Z))^2 FILE_PATH $1.dat

# let us compute the sum of the square deviations, there
# are some different ways of doing this indegenously in wasora,
# but we choose to use the good old awk:
SHELL "awk '\{s += \$6\}END\{print \"$1’s squared deviation is\", s\}' $1.dat"

# the following steps are not important, see the full source tree
# to find out exactly what compare.gp contains

# build an appropriate gnuplot script from a template
SHELL "sed s/xxxxx/$1/ compare.gp > $1.gp"
# and make figures!
SHELL "gnuplot $1.gp"
$ wasora compare.was the14
the14’s squared deviation is 31.9224
$ gnuplot theoretical.gp
$ wasora compare.was seg65
seg65’s squared deviation is 279.774
$ wasora compare.was mey67
mey67’s squared deviation is 263.929
$ wasora compare.was coh71
coh71’s squared deviation is 938.892
$ wasora compare.was pea86
pea86’s squared deviation is 88.4406
$ wasora compare.was the05
the05’s squared deviation is 97.7507
$ 

The fit obtained by wasora gives the following results when \(W(A,Z)\) is evaluated for the 2011 nuclides contained in the experimental set is practically indistinguishable from the actual experimental data, which are shown again:

FIGURE: [028-mass/theoretical2d] Theoretical binding energy for 2011 nuclides given by the semi-empirical mass formula fitted with wasora.

FIGURE: [028-mass/experimental2d] Experimental binding energy for 2011 nuclides. Data taken from the “2003 Atomic Mass Evaluation”: “http://www.nndc.bnl.gov/masses/”

FIGURE: [028-mass/theoretical3d] Surface view of theoretical binding energy for 2011 nuclides given by the semi-empirical mass formula fitted with wasora.

FIGURE: [028-mass/experimental3d] Surface view of experimental binding energy for 2011 nuclides. Data taken from the “2003 Atomic Mass Evaluation”: “http://www.nndc.bnl.gov/masses/”

To effectively compare the fit to the data, its difference is to be analyzed. The figures that follow show the difference between the experimental data and the different fits with a color scale. The greener the plot, the smaller the difference:

FIGURE: [028-mass/seg652d] Difference as computed by Segre’s coefficients

FIGURE: [028-mass/coh712d] Difference as computed by Cohen’s coefficients (probably uses another form of the formula)

FIGURE: [028-mass/mey672d] Difference as computed by Meyerhoff’s coefficients (probably uses another form of the formula)

FIGURE: [028-mass/pea862d] Difference as computed by Pearson’s coefficients

FIGURE: [028-mass/the052d] Difference as computed by me back in 2005

FIGURE: [028-mass/the142d] Difference as computed by wasora

Same information from a different point of view:

FIGURE: [028-mass/seg653d] Difference as computed by Segre’s coefficients

FIGURE: [028-mass/coh713d] Difference as computed by Cohen’s coefficients (probably uses another form of the formula)

FIGURE: [028-mass/mey673d] Difference as computed by Meyerhoff’s coefficients (probably uses another form of the formula)

FIGURE: [028-mass/pea863d] Difference as computed by Pearson’s coefficients

FIGURE: [028-mass/the053d] Difference as computed by me back in 2005

FIGURE: [028-mass/the143d] Difference as computed by wasora

Future versions of “wasora”: “http://www.talador.com.ar/jeremy/wasora” may allow the inclusion of uncertainties or weighting factors associated to the experimental data, so one may pay more attention to either stable nuclei or to those nuclei whose data is more reliable.

Of course, the coefficients reported in nuclear physics textbook may refer to other versions of Weiszacker’s formula, may be fitted against another set of experimental data or whatever. Please do not understand these results as a claim that Dr.~Segre was wrong and I am not (in fact he was not!). Take this report as an example of how “wasora”: “http://www.talador.com.ar/jeremy/wasora” can be used to fit non-linear multidimensional data. Actually, it is less than that. This report is a new chapter of my own scientific history.