Merge pull request #16 from samwaseda/correct_joss_paper

Correct joss paper

CONTRIBUTING.md

@@ -0,0 +1,14 @@
+# Contributing to mamonca
+
+Welcome to mamonca! We appreciate your interest in contributing.
+
+## How to Contribute
+
+1. Fork the repository to your GitHub account.
+2. Clone the forked repository to your local machine:
+
+   ```bash
+   git clone https://github.com/samwaseda/mamonca.git
+   ```
+
+Feel free to submit any changes.

README.md

@@ -2,6 +2,14 @@
 
 This code allows you to launch Metropolis Monte Carlo simulations via Heisenberg Landau models (with various polynomial degrees) from a jupyter notebook.
 
+## Model
+
+`mamonca` is based on the Heisenberg Landau model of the format:
+
+$$\mathcal H = -\frac{1}{2}\sum_{ij,\kappa}J_{ij,\kappa}m_i^{2\kappa+1}m_j^{2\kappa+1} + \sum_{i,\kappa} A_{i,\kappa} m_i^{2\kappa}$$
+
+where $i$ and $j$ go over all atoms and $\kappa$ is the exponent ($\kappa=1$ and $A_{i,\kappa}=0$ for all $i$ and $\kappa$ for the classical Heisenberg model). The evaluation takes place either via Metropolis Monte Carlo (MC) method or spindynamics (SD). MC has the advantage of converging very fast, while SD also delivers the kinetics.
+
 ## How to compile
 
 `mamonca` can be installed directly from conda:
@@ -18,18 +26,19 @@ python setup.py build_ext --user
 
 ## First steps:
 
-In the following simple (but complete) example, we create a bcc Fe system using [pyiron](http://github.com/pyiron/pyiron) and launch a Metropolis Monte Carlo simulation with a Heisenberg coefficient `J=0.1` (eV) for the first nearest neighbor pairs:
+In the following simple (but complete) example, we create a bcc Fe system using [pyiron](http://github.com/pyiron/pyiron) (install via `conda install pyiron`) and launch a Metropolis Monte Carlo simulation with a Heisenberg coefficient `J=0.1` (eV) for the first nearest neighbor pairs:
 
 ```python
 from pyiron_atomistics import Project
 from mamonca import MC
 
-structure = Project('.').create.structure.bulk(
+basis = Project('.').create.structure.bulk(
     name='Fe',
     cubic=True
 )
 
-structure.set_repeat(10)
+# Repeat the structure 10 times in each direction
+structure = basis.repeat(10)
 J = 0.1 # eV
 neighbors = structure.get_neighbors()
 first_shell_tensor = neighbors.get_shell_matrix()[0]
@@ -40,6 +49,13 @@ mc.set_heisenberg_coeff(J * first_shell_tensor)
 mc.run(temperature=300, number_of_iterations=1000)
 ```
 
+More complete list of examples can be found in `notebooks/first_steps.ipynb`
+
 ## How to set inputs and get outputs
 
-As a rule of thumb, you can set all input parameters via functions starting with `set_`. Similarly, output values can be obtained via functions whose names start with `get_`. Take a look at the list of auto-complete and see their docstrings
+As a rule of thumb, you can set all input parameters via functions starting with `set_`. Similarly, output values can be obtained via functions whose names start with `get_`. Most notably, you can get all basic output values via `get_output()` in a dictionary. Otherwise, take a look at the list of auto-complete and see their docstrings
+
+## Notes
+
+- Currently only Linux installation is supported
+- You can run tests located in the `tests` folder

mamonca/cMC.cpp

@@ -527,26 +527,30 @@ void cMC::run_mc(double kBT){
         id_rand = selectable_id.at(rand()%selectable_id.size());
         atom[id_rand].propose_new_state();
         double dEE = atom[id_rand].dE();
+        // Append metadynamics energy if defined
         if (meta.initialized)
             dEE += meta.get_biased_energy(
                 m_norm(magnetization+atom[id_rand].delta_m()/(double)n_tot),
                 sqrt((magnetization*magnetization).sum()));
+        // Suggest a new state
         update_magnetization(id_rand);
         if(thermodynamic_integration())
             dEE = (1-lambda)*dEE+lambda*atom[id_rand].dE(1);
-        if(metropolis(kBT, dEE))
+        // Metropolis step
+        if(metropolis(kBT, dEE)) // Accepted
         {
-            acc++;
+            acc++; // Acceptance ratio
             dEE_tot.at(0) += atom[id_rand].dE();
             if(thermodynamic_integration())
                 dEE_tot.at(1) += atom[id_rand].dE(1);
         }
-        else
+        else /* Rejected */
         {
             update_magnetization(id_rand, true);
             atom[id_rand].revoke();
         }
     }
+    // Energy consistency only for debugging
     for(int i=0; debug_mode && i<2; i++)
     {
         EE_tot[i] += get_energy(i);

mamonca/cMC.h

@@ -42,8 +42,13 @@ struct Magnitude;;
 
 class Atom{
     private:
+        /* Some variables have two values for the two potentials of the      */
+        /* thermodynamic integration. When thermodynamic integration is not  */
+        /* activated, only the first value is used.                          */
         double mabs, mabs_tmp, E_current[2], dE_current[2], dm, dphi, mmax;
         valarray<double> gradient;
+        /* Heisenberg coefficients (i.e. pairwise interactions) and Landau   */
+        /* coefficients (i.e. on-site coefficients)                          */
         vector<double> heisen_coeff[2], landau_coeff[2];
         vector<Magnitude*> landau_func[2];
         vector<Product*> heisen_func[2];
@@ -78,8 +83,13 @@ class Atom{
         void clear_heisenberg_coeff(int);
         void activate_debug();
         void propose_new_state();
+        /* This is the function that defines the maximum magnitude of the    */
+        /* magnetic moment length and the angle change at each step. These   */
+        /* values should be adjusted to have a good acceptance ratio (around */
+        /* 0.33?)                                                            */
         void rescale_magnitude(double, double);
         void set_magnitude(double, double, bool flip_in=true);
+        // Consistency check only in the debugging mode
         void check_consistency();
 };
 

notebooks/first_steps.ipynb

@@ -121,7 +121,7 @@
    "id": "38c184fb",
    "metadata": {},
    "source": [
-    "Here, as we chose the parameter randomly, the Curie temperature is somewhat higher than what it is supposed to be. In addition, since the system is relatively small, the total magnetization after the transition does not completely vanish."
+    "Here, as we chose the Heisenberg parameter randomly, the Curie temperature is somewhat higher than the experimental value of 1,043 K. In addition, since the system is relatively small, the total magnetization after the transition does not completely vanish. The residual magnetization values above the Curie temperature can therefore vary each time you run the simulation. The unit if the magnetization is given by the Heisenberg parameter."
    ]
   },
   {
@@ -132,6 +132,14 @@
     "## Thermodynamic integration"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "ef3e00f3-6e53-44d7-ab4e-3ebd60133215",
+   "metadata": {},
+   "source": [
+    "Let's now take a look at the energy difference of iron in the two most prominent phases: Face-Centered-Cubic (fcc) and Body-Centered-Cubic (bcc). The energy difference can be calculated using [Thermodynamic Integration](https://en.wikipedia.org/wiki/Thermodynamic_integration), where we vary the $\\lambda$ parameter (cf. [Wikipedia](https://en.wikipedia.org/wiki/Thermodynamic_integration#Derivation)) and integrate over the measured potential energies."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 37,
@@ -228,20 +236,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 53,
+   "execution_count": 18,
    "id": "b70a0bb0",
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "The free energy difference between bcc and fcc at 300 K is 0.956543594612767\n"
+      "The free energy difference between bcc and fcc at 300 K is 0.9558411678020047 eV\n"
      ]
     }
    ],
    "source": [
-    "print(\"The free energy difference between bcc and fcc at 300 K is\", E_diff.sum() * np.diff(ti_lambda).mean())"
+    "print(\"The free energy difference between bcc and fcc at 300 K is\", E_diff.sum() * np.diff(ti_lambda).mean(), \"eV\")"
    ]
   },
   {
@@ -252,6 +260,14 @@
     "## Metadynamics"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "a992c792-8fe8-41be-ba54-cc8f726cee5f",
+   "metadata": {},
+   "source": [
+    "[Metadynamics](https://en.wikipedia.org/wiki/Metadynamics) is a simulation method which allows for the free energy distribution along a collective variable defined by the user. In mamonca, you can use the total magnetization as the collective variable, meaning you can obtain the free energy distribution along the total magnetization. For this, you only need to call `set_metadynamics` before `run`, where you must also specify the maximum magnetization value. In order to obtain a reliable result, you should also test the robustness by varying `energy_increment` and `length_scale` (cf. [Wikipedia](https://en.wikipedia.org/wiki/Metadynamics#Algorithm), where `energy_increment` is called $\\omega$ and `length_scale` is called $\\sigma$)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 54,
@@ -316,10 +332,18 @@
     "plt.plot(meta[\"magnetization\"], meta[\"free_energy\"])"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "b8bf7959-f89f-43a1-ba9f-ca4e30678c2d",
+   "metadata": {},
+   "source": [
+    "The free energy minimum shows the most stable state. In this case it is around 0.5. Metadynamics, however, often requires a meticulous sampling, meaning with the standard energy increment of 0.001, you might overshoot and miss the global minimum. For a real measurement, it is recommended to make it as small as possible, meaning you should estimate the amount of time needed to run your calculation with Metadynamics with a small energy increment value and see what would be an acceptable total computation time."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "55339279",
+   "id": "9bf9c044-1a27-45b1-a582-cf72fddc5bb8",
    "metadata": {},
    "outputs": [],
    "source": []
@@ -341,7 +365,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.6"
+   "version": "3.11.8"
   }
  },
  "nbformat": 4,

paper.md

@@ -30,11 +30,11 @@ aas-journal: Astrophysical Journal <- The name of the AAS journal.
 
 # Summary
 
-Magnetic interactions account for a significant portion of free energy in certain materials, ranging from iron, whose ground state would be wrongly predicted without considering magnetic interactions [@friak2001ab], to the magnetocaloric effects of Heusler alloys [@weiss1917phenomene], whose magnetic properties could allow for the development of highly efficient refrgiration systems. In materials science, the Heisenberg model is frequently employed to heuristically compute the potential energy. There are two main methods to make use of the Heisenberg model at finite temperature: one is the Monte Carlo method for an efficient free energy minimization, the other is spin dynamics for the calculation of spin configuration evolution.
+Magnetic interactions account for a significant portion of free energy in certain materials, ranging from relatively simple systems such as iron to complex magnetocaloric effects of Heusler alloys [@weiss1917phenomene]. More specifically, in the case of iron, the ground state would be wrongly predicted without considering magnetic interactions [@friak2001ab]. In Heusler systems, the understanding of magnetic properties could allow for the development of highly efficient refrgiration systems. In materials science, the Heisenberg model is frequently employed to heuristically compute the magnetic part of the potential energy. There are two main methods to make use of the Heisenberg model at finite temperature: one is the Monte Carlo method for an efficient free energy minimization, the other is spin dynamics for the calculation of spin configuration evolution. The Monte Carlo method has the advantage of obtaining the free energy rapidly, while spin dynamics would deliver also the kinetics of the system. `mamonca` allows for the evaluation of the Heisenberg Hamiltonian with extended terms using both Monte Carlo method and spin dynamics.
 
 # Statement of need
 
-`mamonca` is a C++-based python software package for the computation of magnetic interactions in solid materials. All inputs and outputs are given by setters (starting with `set_`) and getters (starting with `get_`), in order for `mamonca` to spare file-reading and writing, in strong contrast to other existing software packages [@kawamura2017quantum; @bauer2011alps; @evans2014atomistic; @hellsvik2011uppsala]. As a result, it has excellent interactivity, as the parameters can be changed on-the-fly, as well as the outputs can be retrieved at any interval chosen by the user. With `mamonca`, the user can analyse any structure that can be defined by other software packages such as Atomic Structure Environment (ASE) [@larsen2017atomic] or pyiron [@janssen2019pyiron], as it takes only the exchange parameters and does not require the knowledge of the structure, which is a strong contrast to existing software packages [@kawamura2017quantum; @bauer2011alps]. `mamonca` has also high flexibility in defining the Hamiltonian, as it allows the user to define not only the classical Heisenberg model, but higher order components including the longitudinal variation, as it has been employed for Fe-Mn systems [@schneider2021ab]. The input parameters for the Hamiltonian can be straightforwardly obtained using a workflow tool such as pyiron, or other calculation software packages such as TB2J [@he2021tb2j]. In addition to the classical Monte Carlo and spin-dynamics, `mamonca` allows also for an addition of Metadynamics [@theodoropoulos2000coarse] and magnetic thermodynamic integration [@frenkel2023understanding], which can deliver the free energy variation. To authors' knowledge, it is the only one code that is able to run Monte Carlo calculations with Metadynamics and magnetic thermodynamic integration.
+`mamonca` is a C++-based python software package for the computation of magnetic interactions in solid materials. All inputs and outputs are given by setters (starting with `set_`) and getters (starting with `get_`), in order for `mamonca` to spare file-reading and writing, in strong contrast to other existing software packages [@kawamura2017quantum; @bauer2011alps; @evans2014atomistic; @hellsvik2011uppsala]. As a result, it has excellent interactivity, as the parameters can be changed on the fly, as well as the outputs can be retrieved at any interval chosen by the user. With `mamonca`, the user can analyse any structure that can be defined by other software packages such as Atomic Structure Environment (ASE) [@larsen2017atomic] or pyiron [@janssen2019pyiron], as `mamonca` takes only the exchange parameters and does not require the knowledge of the structure, which is a strong contrast to existing software packages [@kawamura2017quantum; @bauer2011alps]. `mamonca` has also high flexibility in defining the Hamiltonian, as it allows the user to define not only the classical Heisenberg model, but higher order components including the longitudinal variation, as it has been employed for Fe-Mn systems [@schneider2021ab]. The input parameters for the Hamiltonian can be straightforwardly obtained using a workflow tool such as pyiron, or other calculation software packages such as TB2J [@he2021tb2j]. A typical workflow with pyiron would consist of a general set of physical parameters (chemical element, lattice parameter etc.), which is then evaluated by the software of user's choice. The results can be straightforwardly evaluated to obtain the exchange parameters with the existing tools inside pyiron. Finally, `mamonca` can run to deliver the finite temperature effects of the magnetic part. This means, the user in principle needs only to insert physical parameters to obtain the magnetic finite temperature behaviour they are interested in. In addition to the classical Monte Carlo and spin-dynamics, `mamonca` allows also for an addition of Metadynamics [@theodoropoulos2000coarse] and magnetic thermodynamic integration [@frenkel2023understanding], which can deliver the free energy variation. It is crucial to include these features within the code, as they have to be applied at each step of the simulation and cannot be evaluated in the post-processing. To authors' knowledge, it is the only one code that is able to run Monte Carlo calculations with Metadynamics and magnetic thermodynamic integration.
 
 # Acknowledgements
 

setup.py

@@ -26,5 +26,5 @@
     cmdclass={"build_ext": build_ext},
     ext_modules=cythonize([ext], language_level="3"),
     options={'build': {'build_lib': 'mamonca'}},
-    install_requires=["numpy"],
+    install_requires=["numpy", "cython"],
 )