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

@@ -83,7 +83,7 @@ double TransitionMatrix::totalRate() {
  * @param  dx:   Logarithmic step of the grid (default = 0.005)
  * @retval
  */
-Atom::Atom(int Z, double m, int A, NuclearRadiusModel radius_model, double fc,
+Atom::Atom(int Z, double m, int A, double radius, NuclearRadiusModel radius_model, double fc,
            double dx) {
   // Set the properties
   this->Z = Z;
@@ -112,24 +112,36 @@ Atom::Atom(int Z, double m, int A, NuclearRadiusModel radius_model, double fc,
   // Define radius
   if (A == -1) {
     R = -1;
-  } else {
+  } else if (radius < 0) {
     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) {
     V_coulomb = new CoulombFermi2Potential(Z, R, A);
-  } else {
+  } else if (radius_model == SPHERE) {
     V_coulomb = new CoulombSpherePotential(Z, R);
   }
 
@@ -325,9 +337,9 @@ double Atom::sphereNuclearModel(int Z, int A) {
  * (default = -1)
  * @retval
  */
-DiracAtom::DiracAtom(int Z, double m, int A, NuclearRadiusModel radius_model,
+DiracAtom::DiracAtom(int Z, double m, int A, double R, NuclearRadiusModel radius_model,
                      double fc, double dx, int ideal_minshell)
-  : Atom(Z, m, A, radius_model, fc, dx) {
+  : Atom(Z, m, A, R, radius_model, fc, dx) {
   restE = mu * pow(Physical::c, 2);
   LOG(DEBUG) << "Rest energy = " << restE / Physical::eV << " eV\n";
   idshell = ideal_minshell;
@@ -1035,7 +1047,7 @@ TransitionMatrix DiracAtom::getTransitionProbabilities(int n1, int l1, bool s1,
   return tmat;
 }
 
-DiracIdealAtom::DiracIdealAtom(int Z, double m, int A,
+DiracIdealAtom::DiracIdealAtom(int Z, double m, int A, double R,
                                NuclearRadiusModel radius_model, double fc,
                                double dx)
-  : DiracAtom(Z, m, A, radius_model, fc, dx, 1) {}
+  : DiracAtom(Z, m, A, R, radius_model, fc, dx, 1) {}

lib/atom.hpp

@@ -89,17 +89,17 @@ class Atom {
   double static sphereNuclearModel(int Z, int A);
 
  protected:
- 
+
   // Fundamental properties
   int Z, A;                  // Nuclear charge and mass number
   double m, M, mu;           // Mass of the orbiting particle (e.g. muon, electron), of the nucleus, and effective mass of the system
   double R;                  // Nuclear radius
-  NuclearRadiusModel rmodel; // Nuclear model 
+  NuclearRadiusModel rmodel; // Nuclear model
 
   //Grid
   double rc = 1.0;   // Central radius
-  double dx = 0.005; // Step 
-  
+  double dx = 0.005; // Step
+
   //Potential
   CoulombSpherePotential *V_coulomb;
 
@@ -110,7 +110,7 @@ class Atom {
   EConfPotential V_econf;
 
  public:
-  Atom(int Z = 1, double m = 1, int A = -1,
+  Atom(int Z = 1, double m = 1, int A = -1, double radius = -1,
        NuclearRadiusModel radius_model = POINT, double fc = 1.0,
        double dx = 0.005);
 
@@ -188,7 +188,7 @@ class DiracAtom : public Atom {
   double in_eps = 1e-5;
   int min_n = 1000;
 
-  DiracAtom(int Z = 1, double m = 1, int A = -1,
+  DiracAtom(int Z = 1, double m = 1, int A = -1, double R = -1,
             NuclearRadiusModel radius_model = POINT, double fc = 1.0,
             double dx = 0.005, int ideal_minshell = -1);
 
@@ -221,7 +221,7 @@ class DiracAtom : public Atom {
 // the analytical hydrogen-like solution
 class DiracIdealAtom : public DiracAtom {
  public:
-  DiracIdealAtom(int Z = 1, double m = 1, int A = -1,
+  DiracIdealAtom(int Z = 1, double m = 1, int A = -1, double R = -1,
                  NuclearRadiusModel radius_model = POINT, double fc = 1.0,
                  double dx = 0.005);
 };

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
@@ -49,7 +51,7 @@ MuDiracInputFile::MuDiracInputFile() : InputFile() {
   this->defineIntNode("max_nodes_iter", InputNode<int>(100));    // Max iterations in nodes search
   this->defineIntNode("max_state_iter", InputNode<int>(100));    // Max iterations in state search
   this->defineIntNode("uehling_steps", InputNode<int>(100));     // Uehling correction integration steps
-  this->defineIntNode("xr_print_precision", InputNode<int>(-1)); // Number of digits to print out in values in .xr.out file 
+  this->defineIntNode("xr_print_precision", InputNode<int>(-1)); // Number of digits to print out in values in .xr.out file
   this->defineIntNode("state_print_precision", InputNode<int>(-1)); // Number of digits to print out in values in Dirac state output .{state_name}.out files
   this->defineIntNode("verbosity", InputNode<int>(1));           // Verbosity level (1 to 3)
   this->defineIntNode("output", InputNode<int>(1));              // Output level (1 to 3)
@@ -76,8 +78,11 @@ MuDiracInputFile::MuDiracInputFile() : InputFile() {
 DiracAtom MuDiracInputFile::makeAtom() {
   // Now extract the relevant parameters
   int Z = getElementZ(this->getStringValue("element"));
+  double R = 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, R, nucmodel, 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;
 }