Merge pull request #31 from MichaelHeines/inputRadius

Input radius and skin thickness

docs/source/example.rst

@@ -70,4 +70,30 @@ Note the significant changes - the energies are almost three times smaller than
 
 1. try adding more :literal:`xr_lines`, for example L1-M2 and L1-M3;
 2. try adding a range of lines; this can be written as K1:M5-K1:M5. It will compute all transitions within the given ranges that obey the selection rules to be allowed;
-3. try plotting the spectra in the :literal:`.spec.dat files`, using Gnuplot or importing them in software like Excel or Origin.
+3. try plotting the spectra in the :literal:`.spec.dat files`, using Gnuplot or importing them in software like Excel or Origin.
+
+You can also vary the size and shape of the nucleus using the radius and tFermi keywords. For example, you can evaluate what happens if you increase the radius by 1%, or change the skin thickness parameter in the Fermi model to 2.1 fm (rather than the conventional 2.3 fm). Increasing the radius and decreasing the skin thickness should both reduce the transition energies.
+
+.. code-block:: bash
+
+    element: Au
+    isotope: 197
+    xr_lines: K1-L2,K1-L3
+    write_spec: T
+    nuclear_model: FERMI2
+    uehling_correction: T
+    electronic_config: Au
+    radius: 7.09
+    tFermi: 2.5
+
+Which changes the output to the following
+
+.. code-block:: bash
+
+    # Z = 79, A = 197 amu, m = 206.768 au
+    Line    DeltaE (eV)     W_12 (s^-1)
+    K1-L2   5.54813e+06             1.60217e+18
+    K1-L3   5.71535e+06             1.74712e+18
+    # Z = 79, A = 197 amu, m = 206.768 au
+
+The transition energies have indeed shifted down. Even such a minor change has induced a 45 keV shift.

docs/source/keywords.rst

@@ -40,7 +40,9 @@ Floating point keywords
 ~~~~~~~~~~~~~~~~~~~~~~~~
 These keywords accept a non-integer number. It can be written normally (e.g. 105.3) or in scientific notation (e.g. 1.053E2).
 
-* :literal:`mass`:  mass of the particle in atomic units (1 = mass of the electron). By default it's the mass of the muon, 206.7683.
+* :literal:`mass`: mass of the particle in atomic units (1 = mass of the electron). By default it's the mass of the muon, 206.7683.
+* :literal:`radius`: Solid sphere equivalent radius. (rms radius * sqrt(5/3)). By default it's set to values from Angeli and Marinova (2013).
+* :literal:`tFermi`: Skin thickness for a Fermi model nucleus. By default it's 2.3 fm.
 * :literal:`energy_tol`: absolute tolerance for energy convergence when searching for eigenvalues. Iterations will stop once the energy change is smaller than this number, in atomic units. Default is 1E-7.
 * :literal:`energy_damp`: a damping parameter used in steepest descent energy search to ease convergence. Used to multiply the suggested step :math:`\delta E` and make it smaller. Helps avoiding overshooting; fine-tuning it might help to converge difficult calculations, while making it bigger might make convergence faster in simple ones. Default is 0.5.
 * :literal:`max_dE_ratio`: maximum ratio between energy step, :math:`\delta E`, and current energy :math:`E` in convergence search. If the suggested step exceeds this ratio times the guessed energy, it will be rescaled. This also serves as a measure to avoid overshooting and can be tweaked to get around cases of bad convergence. Default is 0.1.

lib/atom.cpp

@@ -77,14 +77,15 @@ double TransitionMatrix::totalRate() {
  * @param  m: Mass of the orbiting particle (e.g. electron)
  * @param  A: Atomic mass (amus, ignored if -1)
  * @param  radius_model: Which NuclearRadiusModel to use
+ * @param  radius: Solid sphere equivalent radius (fm, ignored if -1)
  * @param  fc:   Central point of the grid (corresponding to i = 0), as a
  * fraction of 1/(Z*mu), the 1s orbital radius for this atom, or of the nuclear
  * radius, depending on which one is bigger (default = 1)
  * @param  dx:   Logarithmic step of the grid (default = 0.005)
  * @retval
  */
-Atom::Atom(int Z, double m, int A, NuclearRadiusModel radius_model, double fc,
-           double dx) {
+Atom::Atom(int Z, double m, int A, NuclearRadiusModel radius_model,
+           double radius, double fc, double dx) {
   // Set the properties
   this->Z = Z;
   this->A = A;
@@ -112,19 +113,31 @@ Atom::Atom(int Z, double m, int A, NuclearRadiusModel radius_model, double fc,
   // Define radius
   if (A == -1) {
     R = -1;
-  } else {
+  } else if (radius == -1) {
     switch (radius_model) {
       case POINT:
         R = -1;
         break;
       case SPHERE:
+        R = sphereNuclearModel(Z, A);
+        LOG(INFO) << "SPHERE: R extracted as " << R << "\n";
+        break;
       case FERMI2:
         R = sphereNuclearModel(Z, A);
+        LOG(INFO) << "FERMI2: R extracted as " << R << "\n";
         break;
       default:
         R = -1;
         break;
     }
+  } else {
+    R = radius * Physical::fm;
+  }
+
+  if (radius_model==POINT) {
+    R = -1;
+  } else {
+    LOG(INFO) << "R = " << radius << " fm\n";
   }
 
   if (radius_model == FERMI2) {
@@ -316,6 +329,7 @@ double Atom::sphereNuclearModel(int Z, int A) {
  * @param  m: Mass of the orbiting particle (e.g. electron)
  * @param  A: Atomic mass (amus, ignored if -1)
  * @param  radius_model: Which NuclearRadiusModel to use
+ * @param  radius: Solid sphere equivalent radius (fm, ignored if -1)
  * @param  fc:   Central point of the grid (corresponding to i = 0), as a
  * fraction of 1/(Z*mu), the 1s orbital radius for this atom, or of the nuclear
  * radius, depending on which one is bigger (default = 1)
@@ -326,8 +340,8 @@ double Atom::sphereNuclearModel(int Z, int A) {
  * @retval
  */
 DiracAtom::DiracAtom(int Z, double m, int A, NuclearRadiusModel radius_model,
-                     double fc, double dx, int ideal_minshell)
-  : Atom(Z, m, A, radius_model, fc, dx) {
+                     double radius, double fc, double dx, int ideal_minshell)
+  : Atom(Z, m, A, radius_model, radius, fc, dx) {
   restE = mu * pow(Physical::c, 2);
   LOG(DEBUG) << "Rest energy = " << restE / Physical::eV << " eV\n";
   idshell = ideal_minshell;
@@ -1035,7 +1049,6 @@ TransitionMatrix DiracAtom::getTransitionProbabilities(int n1, int l1, bool s1,
   return tmat;
 }
 
-DiracIdealAtom::DiracIdealAtom(int Z, double m, int A,
-                               NuclearRadiusModel radius_model, double fc,
-                               double dx)
-  : DiracAtom(Z, m, A, radius_model, fc, dx, 1) {}
+DiracIdealAtom::DiracIdealAtom(int Z, double m, int A, NuclearRadiusModel radius_model,
+                               double radius, double fc,double dx)
+  : DiracAtom(Z, m, A, radius_model, radius, fc, dx, 1) {}

lib/atom.hpp

@@ -110,9 +110,8 @@ class Atom {
   EConfPotential V_econf;
 
  public:
-  Atom(int Z = 1, double m = 1, int A = -1,
-       NuclearRadiusModel radius_model = POINT, double fc = 1.0,
-       double dx = 0.005);
+  Atom(int Z = 1, double m = 1, int A = -1, NuclearRadiusModel radius_model = POINT,
+       double radius = -1, double fc = 1.0,double dx = 0.005);
 
   // Basic getters
   double getZ() {
@@ -188,9 +187,8 @@ class DiracAtom : public Atom {
   double in_eps = 1e-5;
   int min_n = 1000;
 
-  DiracAtom(int Z = 1, double m = 1, int A = -1,
-            NuclearRadiusModel radius_model = POINT, double fc = 1.0,
-            double dx = 0.005, int ideal_minshell = -1);
+  DiracAtom(int Z = 1, double m = 1, int A = -1, NuclearRadiusModel radius_model = POINT,
+            double radius = -1, double fc = 1.0, double dx = 0.005, int ideal_minshell = -1);
 
   double getRestE() {
     return restE;
@@ -221,9 +219,8 @@ class DiracAtom : public Atom {
 // the analytical hydrogen-like solution
 class DiracIdealAtom : public DiracAtom {
  public:
-  DiracIdealAtom(int Z = 1, double m = 1, int A = -1,
-                 NuclearRadiusModel radius_model = POINT, double fc = 1.0,
-                 double dx = 0.005);
+  DiracIdealAtom(int Z = 1, double m = 1, int A = -1, NuclearRadiusModel radius_model = POINT,
+                 double radius = -1, double fc = 1.0, double dx = 0.005);
 };
 
 #endif

lib/config.cpp

@@ -28,6 +28,8 @@ MuDiracInputFile::MuDiracInputFile() : InputFile() {
 
   // Double keywords
   this->defineDoubleNode("mass", InputNode<double>(Physical::m_mu));      // Mass of orbiting particle (default: muon mass)
+  this->defineDoubleNode("radius", InputNode<double>(-1));                // Solid sphere equivalent radius
+  this->defineDoubleNode("tFermi", InputNode<double>(-1));                // Skin thickness for Fermi model
   this->defineDoubleNode("energy_tol", InputNode<double>(1e-7));          // Tolerance for electronic convergence
   this->defineDoubleNode("energy_damp", InputNode<double>(0.5));          // "Damping" used in steepest descent energy search
   this->defineDoubleNode("max_dE_ratio", InputNode<double>(0.1));         // Maximum |dE|/E ratio in energy search
@@ -76,8 +78,11 @@ MuDiracInputFile::MuDiracInputFile() : InputFile() {
 DiracAtom MuDiracInputFile::makeAtom() {
   // Now extract the relevant parameters
   int Z = getElementZ(this->getStringValue("element"));
+  double radius = this->getDoubleValue("radius");
+  double t = this->getDoubleValue("tFermi");
   double m = this->getDoubleValue("mass");
   int A = this->getIntValue("isotope");
+
   if (A == -1) {
     A = getElementMainIsotope(Z);
   }
@@ -98,7 +103,7 @@ DiracAtom MuDiracInputFile::makeAtom() {
 
   // Prepare the DiracAtom
   DiracAtom da;
-  da = DiracAtom(Z, m, A, nucmodel, fc, dx, idshell);
+  da = DiracAtom(Z, m, A, nucmodel, radius, fc, dx, idshell);
   da.Etol = this->getDoubleValue("energy_tol");
   da.Edamp = this->getDoubleValue("energy_damp");
   da.max_dE_ratio = this->getDoubleValue("max_dE_ratio");
@@ -122,5 +127,10 @@ DiracAtom MuDiracInputFile::makeAtom() {
                          this->getDoubleValue("econf_rout_min"));
   }
 
+  if (t != -1) {
+    da.setFermi2(t * Physical::fm);
+    LOG(INFO) << "t = " << t << "\n";
+  }
+
   return da;
 }