Add start of MC tutorials

docs/tutorials.rst

@@ -9,3 +9,4 @@ Tutorials
    tutorials/ase
    tutorials/diffpy
    tutorials/neighbors
+   tutorials/mc

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

docs/tutorials/mc.rst

@@ -0,0 +1,9 @@
+Monte Carlo simulations
+=======================
+
+.. toctree::
+   :maxdepth: 2
+
+   ising
+   size_effect
+

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

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
 
 """