Merge branch 'master' into fast_w_update

Author: mesonepigreco

.github/workflows/python-testsuite.yml

@@ -49,5 +49,6 @@ jobs:
       run: |
         export OMP_NUM_THREADS=1
         cd tests
+        rm -rf __pycache__
         # Test excluding very long running tests
         pytest -v -m "not release"

Examples/ThermodynamicsOfGold/Au.cif

@@ -0,0 +1,27 @@
+# generated using pymatgen
+data_Au
+_symmetry_space_group_name_H-M   'P 1'
+_cell_length_a   2.94954603
+_cell_length_b   2.94954603
+_cell_length_c   2.94954603
+_cell_angle_alpha   60.00000000
+_cell_angle_beta   60.00000000
+_cell_angle_gamma   60.00000000
+_symmetry_Int_Tables_number   1
+_chemical_formula_structural   Au
+_chemical_formula_sum   Au1
+_cell_volume   18.14473112
+_cell_formula_units_Z   1
+loop_
+ _symmetry_equiv_pos_site_id
+ _symmetry_equiv_pos_as_xyz
+  1  'x, y, z'
+loop_
+ _atom_site_type_symbol
+ _atom_site_label
+ _atom_site_symmetry_multiplicity
+ _atom_site_fract_x
+ _atom_site_fract_y
+ _atom_site_fract_z
+ _atom_site_occupancy
+  Au  Au0  1  0.00000000  0.00000000  0.00000000  1

Examples/ThermodynamicsOfGold/get_gold_free_energy.py

@@ -0,0 +1,99 @@
+# Import the sscha code
+import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax, sscha.Utilities
+
+# Import the cellconstructor library to manage phonons
+import cellconstructor as CC, cellconstructor.Phonons
+import cellconstructor.Structure, cellconstructor.calculators
+
+# Import the force field of Gold
+import ase, ase.calculators
+from ase.calculators.emt import EMT
+
+# Import numerical and general pourpouse libraries
+import numpy as np, matplotlib.pyplot as plt
+import sys, os
+
+
+
+"""
+Here we load the primitive cell of Gold from a cif file.
+And we use CellConstructor to compute phonons from finite differences.
+The phonons are computed on a q-mesh 4x4x4
+"""
+
+gold_structure = CC.Structure.Structure()
+gold_structure.read_generic_file("Au.cif")
+
+# Get the force field for gold
+calculator = EMT()
+
+# Relax the gold structure (useless since for symmetries it is already relaxed)
+relax = CC.calculators.Relax(gold_structure, calculator)
+gold_structure_relaxed = relax.static_relax()
+
+# Compute the harmonic phonons
+# NOTE: if the code is run with mpirun, the calculation goes in parallel
+gold_harmonic_dyn = CC.Phonons.compute_phonons_finite_displacements(gold_structure_relaxed, calculator, supercell = (4,4,4))
+
+# Impose the symmetries and 
+# save the dynamical matrix in the quantum espresso format
+gold_harmonic_dyn.Symmetrize()
+gold_harmonic_dyn.save_qe("harmonic_dyn")
+
+
+# If the dynamical matrix has imaginary frequencies, remove them
+gold_harmonic_dyn.ForcePositiveDefinite()
+
+"""
+gold_harmonic_dyn is ready to start the SSCHA calculation.
+
+Now let us initialize the ensemble, and the calculation at 300 K.
+We will run a NVT calculation, using 100 configurations at each step
+"""
+
+TEMPERATURE = 300
+N_CONFIGS = 50
+MAX_ITERATIONS = 20
+
+# Initialize the random ionic ensemble
+ensemble = sscha.Ensemble.Ensemble(gold_harmonic_dyn, TEMPERATURE)
+
+# Initialize the free energy minimizer
+minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
+minim.set_minimization_step(0.01) 
+
+# Initialize the NVT simulation
+relax = sscha.Relax.SSCHA(minim, calculator, N_configs = N_CONFIGS,
+                          max_pop = MAX_ITERATIONS)
+
+# Define the I/O operations
+# To save info about the free energy minimization after each step
+ioinfo = sscha.Utilities.IOInfo()
+ioinfo.SetupSaving("minim_info")
+relax.setup_custom_functions(custom_function_post = ioinfo.CFP_SaveAll)
+
+
+# Run the NVT simulation (save the stress to compute the pressure)
+relax.relax(get_stress = True)
+
+# If instead you want to run a NPT simulation, use
+# The target pressure is given in GPa.
+#relax.vc_relax(target_press = 0)
+
+# You can also run a mixed simulation (NVT) but with variable lattice parameters
+#relax.vc_relax(fix_volume = True)
+
+# Now we can save the final dynamical matrix
+# And print in stdout the info about the minimization
+relax.minim.finalize()
+relax.minim.dyn.save_qe("sscha_T{}_dyn".format(TEMPERATURE))
+
+
+
+
+
+
+
+
+
+

Examples/ThermodynamicsOfGold/plot_dispersion.py

@@ -0,0 +1,72 @@
+import cellconstructor as CC, cellconstructor.Phonons, cellconstructor.ForceTensor
+import ase, ase.dft.kpoints
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import matplotlib.colors as colors
+
+import sys, os
+
+NQIRR = 13
+#CMAP = "Spectral_r"
+PATH = "GXWXKGL"
+N_POINTS = 1000
+
+SPECIAL_POINTS = {"G": [0,0,0],
+                  "X": [0, .5, .5],
+                  "L": [.5, .5, .5],
+                  "W": [.25, .75, .5],
+                  "K": [3/8., 3/4., 3/8.]}
+
+# Load the harmonic and sscha phonons
+harmonic_dyn = CC.Phonons.Phonons('harmonic_dyn', NQIRR)
+sscha_dyn = CC.Phonons.Phonons('sscha_T300_dyn', NQIRR)
+
+# Get the band path
+qpath, data = CC.Methods.get_bandpath(harmonic_dyn.structure.unit_cell,
+                                      PATH,
+                                      SPECIAL_POINTS,
+                                      N_POINTS)
+xaxis, xticks, xlabels = data # Info to plot correclty the x axis
+
+# Get the phonon dispersion along the path
+harmonic_dispersion = CC.ForceTensor.get_phonons_in_qpath(harmonic_dyn, qpath)
+sscha_dispersion = CC.ForceTensor.get_phonons_in_qpath(sscha_dyn, qpath)
+
+nmodes = harmonic_dyn.structure.N_atoms * 3
+
+# Plot the two dispersions
+plt.figure(dpi = 150)
+ax = plt.gca()
+
+for i in range(nmodes):
+    lbl=None
+    lblsscha = None
+    if i == 0:
+        lbl = 'Harmonic'
+        lblsscha = 'SSCHA'
+        
+    ax.plot(xaxis, harmonic_dispersion[:,i], color = 'k', ls = 'dashed', label = lbl)
+    ax.plot(xaxis, sscha_dispersion[:,i], color = 'r', label = lblsscha)
+
+# Plot vertical lines for each high symmetry points
+for x in xticks:
+    ax.axvline(x, 0, 1, color = "k", lw = 0.4)
+ax.axhline(0, 0, 1, color = 'k', ls = ':', lw = 0.4)
+
+ax.legend()
+    
+# Set the x labels to the high symmetry points
+ax.set_xticks(xticks)
+ax.set_xticklabels(xlabels)
+
+ax.set_xlabel("Q path")
+ax.set_ylabel("Phonons [cm-1]")
+
+plt.tight_layout()
+plt.savefig("dispersion.png")
+plt.show()
+
+    
+

Examples/ThermodynamicsOfGold/plot_thermal_expansion.py

@@ -0,0 +1,83 @@
+import cellconstructor as CC, cellconstructor.Phonons
+import numpy as np
+import matplotlib.pyplot as plt
+
+import scipy, scipy.optimize
+
+import sys, os
+
+"""
+This simple scripts plot the thermal expansion.
+It either directly plots the file produced at the end of thermal_expansion.py
+or it processes the dynamical matrices generated throughout the minimization.
+
+It also fits the V(T) cuve and estimates the volumetric thermal expansion coefficient
+
+alpha = dV/dT / V
+
+At 300 K
+"""
+
+
+# Load all the dynamical matrices and compute volume
+DIRECTORY = "thermal_expansion"
+FILE = os.path.join(DIRECTORY, "thermal_expansion.dat")
+
+
+# Check if the file with the thermal expansion data exists
+if not os.path.exists( FILE):
+    # If the simulation is not ended, load the volume from the dynamical matrices
+    
+    # Get the dynamical matrices
+    all_dyn_files = [x for x in os.listdir(DIRECTORY) if "sscha" in x and x.endswith("dyn1")]
+    temperatures = [float(x.split("_")[-2].replace("T", "")) for x in all_dyn_files]
+
+    # Now sort in order of temperature
+    sortmask = np.argsort(temperatures)
+    all_dyn_files = [all_dyn_files[x] for x in sortmask]
+    temperatures = np.sort(temperatures)
+
+    volumes = np.zeros_like(temperatures)
+
+    for i, fname in enumerate(all_dyn_files):
+        # Load the dynamical matrix
+        # The full_name means that we specify the name including the tail 1
+        dyn = CC.Phonons.Phonons(os.path.join(DIRECTORY, fname), full_name = True)
+        volumes[i] = dyn.get_volumes()
+
+else:
+    # Load the data from the final data file
+    temperatures, volumes = np.loadtxt(FILE, unpack = True)
+
+
+# Prepare the figure and plot the V(T) from the sscha data
+plt.figure(dpi = 150)
+plt.scatter(temperatures, volumes, label = "SSCHA data", color = 'r')
+
+# Fit the data with a quadratic curve
+def parabola(x, a, b, c):
+    return a + b*x + c*x**2
+def diff_parab(x, a, b, c):
+    return b + 2*c*x
+
+popt, pcov = scipy.optimize.curve_fit(parabola, temperatures, volumes,
+                                      p0 = [0,0,0])
+
+# Evaluate the volume thermal expansion
+vol_thermal_expansion = diff_parab(300, *popt) / parabola(300, *popt)
+print("Vol thermal expansion: {} x 10^6  K^-1".format(vol_thermal_expansion * 1e6))
+plt.text(0.6, 0.2, r"$\alpha_v = "+"{:.1f}".format(vol_thermal_expansion*1e6)+r"\times 10^6 $ K$^{-1}$",
+         transform = plt.gca().transAxes)
+
+
+# Plot the fit
+_t_ = np.linspace(np.min(temperatures), np.max(temperatures), 1000)
+plt.plot(_t_, parabola(_t_, *popt), label = "Fit", color = 'k', zorder = 0)
+
+# Adjust the plot adding labels, legend, and saving in eps
+plt.xlabel("Temperature [K]")
+plt.ylabel(r"Volume [$\AA^3$]")
+plt.legend()
+plt.tight_layout()
+plt.savefig("thermal_expansion.png")
+plt.show()

Examples/ThermodynamicsOfGold/thermal_expansion.py

@@ -0,0 +1,88 @@
+# Import the sscha code
+import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax, sscha.Utilities
+
+# Import the cellconstructor library to manage phonons
+import cellconstructor as CC, cellconstructor.Phonons
+import cellconstructor.Structure, cellconstructor.calculators
+
+# Import the force field of Gold
+import ase, ase.calculators
+from ase.calculators.emt import EMT
+
+# Import numerical and general pourpouse libraries
+import numpy as np, matplotlib.pyplot as plt
+import sys, os
+
+
+"""
+You need first to run the
+get_gold_free_energy.py
+
+Here we use NPT simulation to compute the gold thermal expansion.
+"""
+
+# Define the temperature range (in K)
+T_START = 300
+T_END = 1000
+DT = 50
+
+N_CONFIGS = 50
+MAX_ITERATIONS = 10
+
+# Import the gold force field
+calculator = EMT()
+
+# Import the starting dynamical matrix (final result of get_gold_free_energy.py)
+dyn = CC.Phonons.Phonons("sscha_T300_dyn", nqirr = 13)
+
+# Create the directory on which to store the output
+DIRECTORY = "thermal_expansion"
+if not os.path.exists(DIRECTORY):
+    os.makedirs("thermal_expansion")
+
+# We cycle over several temperatures
+t = T_START
+
+
+volumes = []
+temperatures = []
+while t <= T_END:
+    # Change the temperature
+    ensemble = sscha.Ensemble.Ensemble(dyn, t)
+    minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
+    minim.set_minimization_step(0.1)
+
+    relax = sscha.Relax.SSCHA(minim, calculator, N_configs = N_CONFIGS,
+                              max_pop = MAX_ITERATIONS)
+
+    # Setup the I/O
+    ioinfo = sscha.Utilities.IOInfo()
+    ioinfo.SetupSaving( os.path.join(DIRECTORY, "minim_t{}".format(t)))
+    relax.setup_custom_functions( custom_function_post = ioinfo.CFP_SaveAll)
+
+
+    # Run the NPT simulation
+    relax.vc_relax(target_press = 0)
+
+    # Save the volume and temperature
+    volumes.append(relax.minim.dyn.structure.get_volume())
+    temperatures.append(t)
+
+    # Start the next simulation from the results
+    relax.minim.dyn.save_qe( os.path.join(DIRECTORY, "sscha_T{}_dyn".format(t)))
+    dyn = relax.minim.dyn
+    relax.minim.finalize()
+    
+    # Update the temperature
+    t += DT
+
+# Save the thermal expansion
+np.savetxt(os.path.join(DIRECTORY, "thermal_expansion.dat"),
+           np.transpose([temperatures, volumes]),
+           header = "Temperature [K]; Volume [A^3]")
+
+    
+
+
+
+