#!/usr/bin/env python # coding: utf-8 # # Quantum relaxation of H3S # Sulfur hydrate is one of the most interesting system of physics, it is one of the highest Tc ever measured (203 K). # Thanks to the very light mass of hydrogen atoms, its physical properties are dominated by anharmonicity and quantum fluctuations, making it the on of the best systems to be tackled with self-consistent harmonic approximation. # # In this simple tutorial we will learn how to run an automatic quantum relaxation. # We will follow the lines of these works: [Errea et al. Nature volume 532, pages81–84(2016)](https://www.nature.com/articles/nature17175) and [Bianco et al Phys. Rev. B 97, 214101](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.97.214101). # # I provide the H3S.scf file, that contains the lattice and the atomic position of the high symmetry phase of H3S. # The harmonic dynamical matrix are given as computed by quantum espresso, however, they can be recomputed live by running the get_phonons.py script (needs ASE and quantum-espresso phonon package). # # You will see that also this case has an instability at gamma, with a phonon mode that is imaginary. # # # How to automatize the relaxation # In the first PbTe tutorial we saw how to prepare the SSCHA calculation from scratch and how to manipulate ensembles and the # minimization. # However, all the procedure could be a bit clumsy. In this section I will show how to automatize everything. # We will complete the PbTe simulation to get some usefull results. # # The automatic relaxation is performed through the Relax module, that implements two classes, one for the static cell relaxation and the other for the variable cell relaxation # In[1]: #get_ipython().run_line_magic('matplotlib', 'notebook') from __future__ import print_function from __future__ import division import sys,os # In[2]: # Import ASE to setup the automatic calculator import ase from ase.calculators.espresso import Espresso # Lets import cellconstructor import cellconstructor as CC import cellconstructor.Phonons # Import the sscha import sscha, sscha.Ensemble, sscha.Relax, sscha.SchaMinimizer from numpy import * import numpy as np import matplotlib.pyplot as plt ## parallel qe using srun (for slurm system) label = 'espresso' input_file = label+'.pwi' output_file = label+'.pwo' no_cpus = 36*4 npool = 36 pw_loc = '/home/adenchfield/packages/q-e-qe-6.5/bin/pw.x' #os.environ['ASE_ESPRESSO_COMMAND'] = ' {} < {} > {}'.format(no_cpus, pw_loc, npool, input_file, output_file) os.environ['ASE_ESPRESSO_COMMAND'] = 'mpirun -np {} {} -npool {} < {} > {}'.format(no_cpus, pw_loc, npool, input_file, output_file) #plt.ion() #Setup interactive plot mode # To make everything clear, we will start again from the harmonic dynamical matrix for our structure: # # matdyn # # We need to impose the sum rule and to force it to be positive definite before starting a SSCHA calculation # In[3]: pseudos = {"H":"H.pbe-rrkjus_psl.1.0.0.UPF", "S" : "S.pbe-nl-rrkjus_psl.1.0.0.UPF"} input_params = {"ecutwfc" : 35, "ecutrho" : 350, "occupations" : "smearing", # "input_dft" : "blyp", "mixing_beta" : 0.2, "conv_thr" : 1e-14, "degauss" : 0.02, "smearing" : "mv", "press" : 2100, "cell_dofree" : "ibrav", "pseudo_dir" : ".", "tstress" : True, "tprnfor" : True} espresso_calc = Espresso(input_data = input_params, calcstress=True, pseudopotentials = pseudos, kpts=(5, 5, 5), koffset=(1,1,1)) # , kspacing = k_spacing, koffset=(1,1,1) # In[4]: # Let us load the starting dynamical matrix dyn = CC.Phonons.Phonons("harmonic_dyn", nqirr = 4) # Apply the sum rule and delete the imaginary modes dyn.Symmetrize() dyn.ForcePositiveDefinite() # Generate the ensemble and the minimizer objects NUM_ENSEMBLES = 100 # 1000 by default ensemble = sscha.Ensemble.Ensemble(dyn, NUM_ENSEMBLES, supercell = dyn.GetSupercell()) minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble) # We setup all the minimization parameters minimizer.min_step_dyn = 0.001 minimizer.kong_liu_ratio = 0.5 # You can automatize also the cluster calculation by setting up a cluster object # like we did in the previous tutorials. # If you keep it as None (as done in the following cell) the calculation will be runned locally. # Remember, if you use the cluster, you need to copy first the pseudopotential in the working directory of the cluster. # In[5]: print("Here") # Here we prepare a cluster # Here we configure the cluster object MARCONI import sscha.Cluster # In[6]: input_params = {"ecutwfc" : 35, "ecutrho" : 350, "occupations" : "smearing", # "input_dft" : "blyp", "mixing_beta" : 0.2, "conv_thr" : 1e-14, "degauss" : 0.02, "smearing" : "mv", "press" : 2100, "cell_dofree" : "ibrav", "pseudo_dir" : ".", "tstress" : True, "tprnfor" : True} input_params["calculation"] = "scf" # Setup the simple DFT calculation espresso_calc = Espresso(input_data = input_params, calcstress=True, pseudopotentials = pseudos, kpts=(5, 5, 5), koffset=(1,1,1)) # , kspacing = k_spacing # , koffset = (1,1,1) # We prepare the Automatic relaxation relax = sscha.Relax.SSCHA(minimizer, ase_calculator = espresso_calc, N_configs = 200, max_pop = 6, save_ensemble = True, cluster = None) import spglib print ("The original spacegroup is:", spglib.get_spacegroup(dyn.structure.get_ase_atoms(), 0.05)) # we define a function that prints the space group during the optimization space_groups = [] def print_spacegroup(minim): spgroup = spglib.get_spacegroup(minim.dyn.structure.get_ase_atoms(), 0.05) space_groups.append(spgroup) # We can save them in the output at each minimization step f = open("space_group.dat", "w") f.writelines(["{}) {}\n".format(i+1, x) for i,x in enumerate(space_groups)]) f.close() relax.setup_custom_functions(custom_function_post = print_spacegroup) print("Here 2") # Lets see the parameters: # # minimizer : it is the SSCHA_Minimizer, containing the settings of each minimization. # N_configs : the number of configurations to generate at each run (we call them populations) # max_pop : The maximum number of populations after wich the calculation is stopped even if not converged. # save_ensemble : If True, after each energy and force calculation, the ensemble will be saved. # cluster : The cluster object to be used for submitting the configurations. # # # If no cluster is provided (like in this case) the calculation is performed locally. # The save_ensemble keyword will dave the ensemble inside the directory specified in the running command # ```python # relax.relax(ensemble_loc = "directory_of_the_ensemble") # ``` # # The calculation can be runned just by calling # ```python # relax.relax() # ``` # And the code will proceed with the automatic SSCHA minimization. # However, before doing so, we will setup a custom function. # These are usefull to manipulate a bit the minimization, printing or saving extra info # during the minimization. # In particular, we will setup a function to be called after each minimization step of each population, that # will save the current frequencies of the auxiliary dynamical matrix in an array, so that we will be able # to plot the whole frequency evolution after the minimization. # # The following cell will start the computation. Consider that it may take long time. # Running espresso in parallel with 4 processrs on a Intel(R) Core(TM) i7-4790K CPU at 4.00GHz the single ab-initio run takes a bit more than 1 minute. In this example we are running 20 configurations per population and up to 6 populations. # So the overall time is about 2 hours and half. # # For hevier relaxations it is strongly suggested to use the cluster module, or to directly run the SSCHA from a more powerful computer. # # Thansk to the line inside our custom function # ```python # np.savetxt("all_frequencies.dat", all_frequencies) # ``` # we will save at each minimization step the evolution of the frequencies. # So we can load and plot them even during the minimization, to get a clue on how frequencies of the SSCHA dynamical matrix are evolving. # # **NOTE**: These frequencies are not the physical frequencies of the phonon quasiparticles, but just the frequencies of the auxiliary dynamical matrix. To obtain the real physical frequencies, look at the StructuralInstability and the Dynamical spectral function tutorials. # In[ ]: # We reset the frequency array all_frequencies = [] # We define a function that will be called by the code after each minimization step, # passing to us the minimizer at that point # We will store the frequencies, to plot after the minimization their evolution. def add_current_frequencies(minimizer): # Get the frequencies w, p = minimizer.dyn.DiagonalizeSupercell() all_frequencies.append(w) # In this way the file will be updated at each step np.savetxt("all_frequencies.dat", all_frequencies) # We add this function to the relax ojbect relax.setup_custom_functions(custom_function_post = add_current_frequencies) # Now we are ready to # ***************** RUN ***************** relax.vc_relax(ensemble_loc = "data_ensemble_autorelax", target_press=2100, static_bulk_modulus=200) # *************************************** print("Here 3") # In[8]: relax.minim.finalize() # We can save the frequency list for furture analisys # NOTE: This will overwrite the one given in the examples! np.savetxt("all_frequencies.dat", all_frequencies) # # In[9]: # Lets plot the minimization relax.minim.plot_results(save_filename="H3S_minimizeroutput", plot=False) print(relax.minim.dyn.structure.unit_cell) # # The Structural phase transition # # We got the minimum of the free energy, now we can compute the free energy curvature around this high symmetry structure, to discover if there are some imaginary phonons. # In order to do this, we need the the free energy hessian: # # $$ # \frac{\partial^2 F}{\partial R_a \partial R_b} \approx \Phi_{ab} + \sum_{pqrs}\stackrel{(3)}{\Phi}_{apq} \Lambda_{pqrs} \stackrel{(3)}{\Phi}_{rsb} # $$ # # This is obtained neglecting 4 phonon scattering (that may be relevant for highly anharmonic systems). # # You can compute this quantity using the get_free_energy_hessian subroutine of the ensemble class. # However, you need a bigger ensemble to get reliable results. We will use the last ensemble, but for production calculations it is advisable to generate a new, bigger, ensemble. # # Here, the flag include_v4 tells if you need to include the V4 # In[12]: print("Right before free_energy_Hessian") # Compute the free energy hessian free_energy_hessian = relax.minim.ensemble.get_free_energy_hessian(include_v4 = False, verbose = True) print("Right after free energy Hessian") # In[14]: # Now we can save the free energy hessian like it was a dynamical matrix free_energy_hessian.save_qe("hessian") # We can print its eigenvalues to see if there are imaginary phonons w, p = free_energy_hessian.DiagonalizeSupercell() print("Hessian eigenvalues") print(w * CC.Units.RY_TO_CM) # No imaginary phonon! The structure is stable at this pressure. # # Remember, to get converged result, you should study this result as a function of the number of configurations and as a function of the supercell size! # We can save the dynamical matrix relax.minim.dyn.save_qe("dyn_pop1_") # Print the frequencies before and after the minimization w_old, p_old = ensemble.dyn_0.DiagonalizeSupercell() # This is the representation of the density matrix used to generate the ensemble w_new, p_new = relax.minim.dyn.DiagonalizeSupercell() # We can now print them print(" Old frequencies | New frequencies") print("\n".join(["{:16.4f} | {:16.4f} cm-1".format(w_old[i] * CC.Units.RY_TO_CM, w_new[i] * CC.Units.RY_TO_CM) for i in range(len(w_old))]))