From a3538193343b5e086c8f4557d1890bac93dfddb8 Mon Sep 17 00:00:00 2001 From: Ross Whitfield Date: Tue, 10 Dec 2019 16:56:23 -0500 Subject: [PATCH] Add start of MC tutorials --- docs/tutorials.rst | 1 + docs/tutorials/ising.rst | 120 +++++++++++++++++++++++++++++++++ docs/tutorials/mc.rst | 9 +++ docs/tutorials/size_effect.rst | 59 ++++++++++++++++ javelin/energies.pyx | 4 +- 5 files changed, 191 insertions(+), 2 deletions(-) create mode 100644 docs/tutorials/ising.rst create mode 100644 docs/tutorials/mc.rst create mode 100644 docs/tutorials/size_effect.rst diff --git a/docs/tutorials.rst b/docs/tutorials.rst index 37e4c82..aab2231 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -9,3 +9,4 @@ Tutorials tutorials/ase tutorials/diffpy tutorials/neighbors + tutorials/mc diff --git a/docs/tutorials/ising.rst b/docs/tutorials/ising.rst new file mode 100644 index 0000000..c3cf930 --- /dev/null +++ b/docs/tutorials/ising.rst @@ -0,0 +1,120 @@ +Ising Model examples +==================== + +Creating chemical short-range order + +Negative correlation +-------------------- + +.. plot:: + :include-source: + + >>> import numpy as np + >>> import xarray as xr + >>> from javelin.structure import Structure + >>> from javelin.energies import IsingEnergy + >>> from javelin.modifier import SwapOccupancy + >>> from javelin.fourier import Fourier + >>> from javelin.mc import MC + >>> n=128 + >>> structure = Structure(symbols=np.random.choice(['Na','Cl'],n**2), positions=[(0., 0., 0.)]*n**2, unitcell=5) + >>> structure.reindex([n,n,1,1]) + >>> e1 = IsingEnergy(11,17,J=0.5) + >>> nl = structure.get_neighbors()[0,-1] + >>> mc = MC() + >>> mc.add_modifier(SwapOccupancy(0)) + >>> mc.temperature = 1 + >>> mc.cycles = 50 + >>> mc.add_target(nl, e1) + >>> out = mc.run(structure) # doctest: +SKIP + >>> f=Fourier() + >>> f.grid.bins = 200, 200, 1 + >>> f.grid.r1 = -3,3 + >>> f.grid.r2 = -3,3 + >>> f.lots = 8,8,1 + >>> f.number_of_lots = 256 + >>> f.average = True + >>> results=f.calc(out) # doctest: +SKIP + >>> # plot + >>> fig, axs = plt.subplots(2, 1, figsize=(6.4,9.6)) # doctest: +SKIP + >>> xr.DataArray.from_series(out.atoms.Z).plot(ax=axs[0]) # doctest: +SKIP + >>> results.plot(ax=axs[1]) # doctest: +SKIP + +Positive correlation +-------------------- + +.. plot:: + :include-source: + + >>> import numpy as np + >>> import xarray as xr + >>> from javelin.structure import Structure + >>> from javelin.energies import IsingEnergy + >>> from javelin.modifier import SwapOccupancy + >>> from javelin.fourier import Fourier + >>> from javelin.mc import MC + >>> + >>> n=100 + >>> structure = Structure(symbols=np.random.choice(['Na','Cl'],n**2), positions=[(0., 0., 0.)]*n**2, unitcell=5) + >>> structure.reindex([n,n,1,1]) + >>> e1 = IsingEnergy(11,17,J=-0.5) + >>> nl = structure.get_neighbors()[0,-1] + >>> mc = MC() + >>> mc.add_modifier(SwapOccupancy(0)) + >>> mc.temperature = 1 + >>> mc.cycles = 50 + >>> mc.add_target(nl, e1) + >>> out = mc.run(structure) # doctest: +SKIP + >>> f=Fourier() + >>> f.grid.bins = 121,121,1 + >>> f.grid.r1 = -3,3 + >>> f.grid.r2 = -3,3 + >>> f.average = True + >>> results=f.calc(out) # doctest: +SKIP + >>> # plot + >>> fig, axs = plt.subplots(2, 1, figsize=(6.4,9.6)) # doctest: +SKIP + >>> xr.DataArray.from_series(out.atoms.Z).plot(ax=axs[0]) # doctest: +SKIP + >>> results.plot(ax=axs[1]) # doctest: +SKIP + + +Getting a desired correlation +----------------------------- + +.. plot:: + :include-source: + + >>> import numpy as np + >>> import xarray as xr + >>> from javelin.structure import Structure + >>> from javelin.energies import IsingEnergy + >>> from javelin.modifier import SwapOccupancy + >>> from javelin.fourier import Fourier + >>> from javelin.mc import MC + >>> + >>> n=100 + >>> structure = Structure(symbols=np.random.choice(['Na','Cl'],n**2), positions=[(0., 0., 0.)]*n**2, unitcell=5) + >>> structure.reindex([n,n,1,1]) + >>> e1 = IsingEnergy(11,17,desired_correlation=0.5) + >>> nl1 = structure.get_neighbors()[0,-1] + >>> e2 = IsingEnergy(11,17,desired_correlation=0) + >>> nl2 = structure.get_neighbors(minD=2.99,maxD=3.01)[0,-1] + >>> e3 = IsingEnergy(11,17,desired_correlation=-0.5) + >>> nl3 = structure.get_neighbors()[1,-2] + >>> mc = MC() + >>> mc.add_modifier(SwapOccupancy(0)) + >>> mc.temperature = 1 + >>> mc.cycles = 50 + >>> mc.add_target(nl1, e1) + >>> mc.add_target(nl2, e2) + >>> mc.add_target(nl3, e3) + >>> out = mc.run(structure) # doctest: +SKIP + >>> f=Fourier() + >>> f.grid.bins = 121,121,1 + >>> f.grid.r1 = -3,3 + >>> f.grid.r2 = -3,3 + >>> f.average = True + >>> results=f.calc(out) # doctest: +SKIP + >>> # plot + >>> fig, axs = plt.subplots(2, 1, figsize=(6.4,9.6)) # doctest: +SKIP + >>> xr.DataArray.from_series(out.atoms.Z).plot(ax=axs[0]) # doctest: +SKIP + >>> results.plot(ax=axs[1]) # doctest: +SKIP diff --git a/docs/tutorials/mc.rst b/docs/tutorials/mc.rst new file mode 100644 index 0000000..0cfb407 --- /dev/null +++ b/docs/tutorials/mc.rst @@ -0,0 +1,9 @@ +Monte Carlo simulations +======================= + +.. toctree:: + :maxdepth: 2 + + ising + size_effect + diff --git a/docs/tutorials/size_effect.rst b/docs/tutorials/size_effect.rst new file mode 100644 index 0000000..2c025fb --- /dev/null +++ b/docs/tutorials/size_effect.rst @@ -0,0 +1,59 @@ +Size-effect +=========== + + +.. plot:: + :include-source: + + >>> import numpy as np + >>> from javelin.structure import Structure + >>> from javelin.energies import LennardJonesEnergy + >>> from javelin.modifier import SetDisplacementNormalXYZ + >>> from javelin.fourier import Fourier + >>> from javelin.mc import MC + >>> n = 128 + >>> + >>> x = np.random.normal(0, 0.01, size=n**2) + >>> y = np.random.normal(0, 0.01, size=n**2) + >>> z = np.zeros(n**2) + >>> + >>> structure = Structure(symbols=np.random.choice(['Na','Cl'],n**2), positions=np.vstack((x,y,z)).T, unitcell=5) + >>> structure.reindex([n,n,1,1]) + >>> + >>> nl = structure.get_neighbors()[0,1,-2,-1] + >>> + >>> e1 = LennardJonesEnergy(1, 1.05, atom_type1=11, atom_type2=11) + >>> e2 = LennardJonesEnergy(1, 1.0, atom_type1=11, atom_type2=17) + >>> e3 = LennardJonesEnergy(1, 0.95, atom_type1=17, atom_type2=17) + >>> + >>> mc = MC() + >>> mc.add_modifier(SetDisplacementNormalXYZ(0, 0, 0.02, 0, 0.02, 0, 0)) + >>> mc.temperature = 0.001 + >>> mc.cycles = 50 + >>> mc.add_target(nl, e1) + >>> mc.add_target(nl, e2) + >>> mc.add_target(nl, e3) + >>> + >>> out = mc.run(structure) # doctest: +SKIP + >>> + >>> f=Fourier() # doctest: +SKIP + >>> f.grid.bins = (201, 201, 1) # doctest: +SKIP + >>> f.grid.r1 = -3, 3 # doctest: +SKIP + >>> f.grid.r2 = -3, 3 # doctest: +SKIP + >>> f.lots = 8, 8, 1 # doctest: +SKIP + >>> f.number_of_lots = 256 # doctest: +SKIP + >>> f.average = True # doctest: +SKIP + >>> + >>> results1=f.calc(out) # doctest: +SKIP + >>> + >>> out.replace_atom(11,42) # doctest: +SKIP + >>> out.replace_atom(17,11) # doctest: +SKIP + >>> out.replace_atom(42,17) # doctest: +SKIP + >>> results2=f.calc(out) # doctest: +SKIP + >>> + >>> out.replace_atom(17,11) # doctest: +SKIP + >>> results3=f.calc(out) # doctest: +SKIP + >>> fig, axs = plt.subplots(3, 1, figsize=(6.4,14.4)) # doctest: +SKIP + >>> results1.plot(ax=axs[0]) # doctest: +SKIP + >>> results2.plot(ax=axs[1]) # doctest: +SKIP + >>> results3.plot(ax=axs[2]) # doctest: +SKIP diff --git a/javelin/energies.pyx b/javelin/energies.pyx index 37d6a82..f134079 100644 --- a/javelin/energies.pyx +++ b/javelin/energies.pyx @@ -20,7 +20,7 @@ with the neighbor vector. a1, x1, y1, z1, a2, x2, y2, z2, neighbor_x, neighbor_y, neighbor_z): - return self.E + return self.E This is slower than using compile classes by about a factor of 10. If you are using IPython or Jupyter notebooks you can use Cython magic to @@ -39,7 +39,7 @@ compile your own energies. You need load the Cython magic first int a1, double x1, double y1, double z1, int a2, double x2, double y2, double z2, Py_ssize_t neighbor_x, Py_ssize_t neighbor_y, Py_ssize_t neighbor_z) except *: - return self.E + return self.E """