From ded91ca0540bc1bfebfd23a361722e04611bf91d Mon Sep 17 00:00:00 2001 From: samwaseda Date: Fri, 23 Feb 2024 09:59:30 +0000 Subject: [PATCH 1/2] Copy files that I created for joss --- CONTRIBUTING.md | 14 ++++++++++++++ README.md | 7 +++++-- mamonca/cMC.cpp | 10 +++++++--- mamonca/cMC.h | 10 ++++++++++ 4 files changed, 36 insertions(+), 5 deletions(-) create mode 100644 CONTRIBUTING.md diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..1191736 --- /dev/null +++ b/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. diff --git a/README.md b/README.md index d5d2d56..f244ec5 100644 --- a/README.md +++ b/README.md @@ -26,7 +26,7 @@ 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 @@ -49,10 +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_`. Most notably, you can get all basic output via `get_output()` in a dictionary. 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 diff --git a/mamonca/cMC.cpp b/mamonca/cMC.cpp index 4e09c2c..888c86f 100644 --- a/mamonca/cMC.cpp +++ b/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); diff --git a/mamonca/cMC.h b/mamonca/cMC.h index 9406fbb..74a5e8a 100644 --- a/mamonca/cMC.h +++ b/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 gradient; + /* Heisenberg coefficients (i.e. pairwise interactions) and Landau */ + /* coefficients (i.e. on-site coefficients) */ vector heisen_coeff[2], landau_coeff[2]; vector landau_func[2]; vector 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(); }; From e10c74964ec38d048ed5e90295b6e82e3bd8bca1 Mon Sep 17 00:00:00 2001 From: samwaseda Date: Tue, 27 Feb 2024 07:48:45 +0000 Subject: [PATCH 2/2] update notebook and setup --- notebooks/first_steps.ipynb | 36 ++++++++++++++++++++++++++++++------ setup.py | 9 +++++++-- 2 files changed, 37 insertions(+), 8 deletions(-) diff --git a/notebooks/first_steps.ipynb b/notebooks/first_steps.ipynb index c3debc4..c05d0b4 100644 --- a/notebooks/first_steps.ipynb +++ b/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 unitless $\\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,7 +236,7 @@ }, { "cell_type": "code", - "execution_count": 53, + "execution_count": 18, "id": "b70a0bb0", "metadata": {}, "outputs": [ @@ -236,12 +244,12 @@ "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, diff --git a/setup.py b/setup.py index 02021f0..6523d35 100644 --- a/setup.py +++ b/setup.py @@ -15,7 +15,7 @@ setup( name='mamonca', - version='0.0.7', + version='0.0.8', description='mamonca - interactive Magnetic Monte Carlo code', long_description=readme, long_description_content_type='text/markdown', @@ -26,5 +26,10 @@ cmdclass={"build_ext": build_ext}, ext_modules=cythonize([ext], language_level="3"), options={'build': {'build_lib': 'mamonca'}}, - install_requires=["numpy", "cython"], + setup_requires=[ + # Setuptools 18.0 properly handles Cython extensions. + 'setuptools>=18.0', + 'cython', + 'numpy', + ], )