From a77844ea1a9bb55e4a40a52f5d3ebf57682c359f Mon Sep 17 00:00:00 2001 From: aichhorn Date: Mon, 30 Mar 2020 15:46:54 +0200 Subject: [PATCH] First draft of the Sr2MgOsO6 w/o SOC doc --- doc/tutorials.rst | 8 + .../images_scripts/Sr2MgOsO6_noSOC.indmftpr | 20 ++ .../images_scripts/Sr2MgOsO6_noSOC.struct | 40 +++ doc/tutorials/sr2mgoso6_nosoc.rst | 268 ++++++++++++++++++ 4 files changed, 336 insertions(+) create mode 100644 doc/tutorials/images_scripts/Sr2MgOsO6_noSOC.indmftpr create mode 100644 doc/tutorials/images_scripts/Sr2MgOsO6_noSOC.struct create mode 100644 doc/tutorials/sr2mgoso6_nosoc.rst diff --git a/doc/tutorials.rst b/doc/tutorials.rst index d02cdeb6..c24133fd 100644 --- a/doc/tutorials.rst +++ b/doc/tutorials.rst @@ -14,6 +14,14 @@ A simple example: SrVO3 tutorials/srvo3 +Basis rotations: Sr2MgOsO6 without SOC +-------------------------------------- + +.. toctree:: + :maxdepth: 2 + + tutorials/sr2mgoso6_nosoc + Full charge self consistency with Wien2k: :math:`\gamma`-Ce ----------------------------------------------------------- diff --git a/doc/tutorials/images_scripts/Sr2MgOsO6_noSOC.indmftpr b/doc/tutorials/images_scripts/Sr2MgOsO6_noSOC.indmftpr new file mode 100644 index 00000000..2f5c5f78 --- /dev/null +++ b/doc/tutorials/images_scripts/Sr2MgOsO6_noSOC.indmftpr @@ -0,0 +1,20 @@ +5 ! Nsort +2 1 1 4 2 ! Mult(Nsort) +3 ! lmax +complex ! choice of angular harmonics +0 0 0 0 ! l included for each sort +0 0 0 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +cubic ! choice of angular harmonics +0 0 2 0 ! l included for each sort +0 0 0 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +0 ! SO flag +complex ! choice of angular harmonics +0 0 0 0 ! l included for each sort +0 0 0 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +complex ! choice of angular harmonics +0 0 0 0 ! l included for each sort +0 0 0 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +complex +0 0 0 0 +0 0 0 0 +-0.088 0.43 ! 0.40 gives warnings, 0.043 gives occ 1.996 diff --git a/doc/tutorials/images_scripts/Sr2MgOsO6_noSOC.struct b/doc/tutorials/images_scripts/Sr2MgOsO6_noSOC.struct new file mode 100644 index 00000000..53439abd --- /dev/null +++ b/doc/tutorials/images_scripts/Sr2MgOsO6_noSOC.struct @@ -0,0 +1,40 @@ +Sr2MgOsO6 +B 5 87 + RELA + 10.507954 10.507954 14.968880 90.000000 90.000000 90.000000 +ATOM -1: X=0.00000000 Y=0.50000000 Z=0.75000000 + MULT= 2 ISPLIT=-2 + -1: X=0.50000000 Y=0.00000000 Z=0.75000000 +Sr 2+ NPT= 781 R0=.000010000 RMT= 2.50000 Z: 38.00000 +LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000 + 0.0000000 1.0000000 0.0000000 + 0.0000000 0.0000000 1.0000000 +ATOM -2: X=0.00000000 Y=0.00000000 Z=0.00000000 + MULT= 1 ISPLIT=-2 +Os 6+ NPT= 781 R0=.000005000 RMT= 1.94 Z: 76.00000 +LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000 + 0.0000000 1.0000000 0.0000000 + 0.0000000 0.0000000 1.0000000 +ATOM -3: X=0.00000000 Y=0.00000000 Z=0.50000000 + MULT= 1 ISPLIT=-2 +Mg 2+ NPT= 781 R0=.000100000 RMT= 1.89 Z: 12.00000 +LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000 + 0.0000000 1.0000000 0.0000000 + 0.0000000 0.0000000 1.0000000 +ATOM -4: X=0.74270000 Y=0.21790000 Z=0.00000000 + MULT= 4 ISPLIT= 8 + -4: X=0.25730000 Y=0.78210000 Z=0.00000000 + -4: X=0.21790000 Y=0.25730000 Z=0.00000000 + -4: X=0.78210000 Y=0.74270000 Z=0.00000000 +O 2- NPT= 781 R0=.000100000 RMT= 1.58 Z: 8.00000 +LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000 + 0.0000000 1.0000000 0.0000000 + 0.0000000 0.0000000 1.0000000 +ATOM -5: X=0.00000000 Y=0.00000000 Z=0.75390000 + MULT= 2 ISPLIT=-2 + -5: X=0.00000000 Y=0.00000000 Z=0.24610000 +O 2- NPT= 781 R0=.000100000 RMT= 1.58 Z: 8.00000 +LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000 + 0.0000000 1.0000000 0.0000000 + 0.0000000 0.0000000 1.0000000 + 0 NUMBER OF SYMMETRY OPERATIONS diff --git a/doc/tutorials/sr2mgoso6_nosoc.rst b/doc/tutorials/sr2mgoso6_nosoc.rst new file mode 100644 index 00000000..36ddfd36 --- /dev/null +++ b/doc/tutorials/sr2mgoso6_nosoc.rst @@ -0,0 +1,268 @@ +.. _Sr2MgOsO6_noSOC: + +Here we will discuss a calculation where off-diagonal matrix elements show up, and will discuss step-by-step how this calculation can be set up. + +The full script for this calculation is also provided here (:download:`dft_dmft_cthyb.py `). + +Note that we do not include spin-orbit coupling here for pedagogical reasons. For the real material it is necessary to include also SOC. + + +DFT (Wien2k) and Wannier orbitals +================================= + +DFT setup +--------- + +First, we do a DFT calculation, using the Wien2k package. As main input file we have to provide the so-called struct file :file:`Sr2MgOs6_noSOC.struct`. We use the following: + +.. literalinclude:: images_scripts/Sr2MgOsO6_noSOC.struct + +The DFT calculation is done as usual, for instance you can use for the initialisation + + + init -b -vxc 5 -numk 2000 + +This is setting up a non-magnetic calculation, using the LDA and 2000 k-points in the full Brillouin zone. As usual, we start the DFT self consistent cycle by the Wien2k script :: + + run + +Wannier orbitals +---------------- + +As a next step, we calculate localised orbitals for the t2g orbitals. We use the same input file for :program:`dmftproj` as it was used in the :ref:`documentation`: + +.. literalinclude:: images_scripts/Sr2MgOsO6_noSOC.indmftpr + +Note that, due to the distortions in the crystal structure, we need to include all five d orbitals in the calculation (line 8 in the input file above). + +To prepare the input data for :program:`dmftproj` we execute lapw2 with the `-almd` option :: + + x lapw2 -almd + +Then :program:`dmftproj` is executed in its default mode (i.e. without spin-polarization or spin-orbit included) :: + + dmftproj + +This program produces the necessary files for the conversion to the hdf5 file structure. This is done using +the python module :class:`Wien2kConverter `. A simple python script that initialises the converter is:: + + from triqs_dft_tools.converters.wien2k_converter import * + Converter = Wien2kConverter(filename = "Sr2MgOsO6_noSOC") + +After initializing the interface module, we can now convert the input +text files to the hdf5 archive by:: + + Converter.convert_dft_input() + +This reads all the data, and stores everything that is necessary for the DMFT calculation in the file :file:`Sr2MgOsO6_noSOC.h5`. + +[CONTINUE HERE] + +The DMFT calculation +==================== + +The DMFT script itself is, except very few details, independent of the DFT package that was used to calculate the local orbitals. +As soon as one has converted everything to the hdf5 format, the following procedure is practially the same. + +Loading modules +--------------- + +First, we load the necessary modules:: + + from triqs_dft_tools.sumk_dft import * + from pytriqs.gf import * + from pytriqs.archive import HDFArchive + from pytriqs.operators.util import * + from triqs_cthyb import * + import pytriqs.utility.mpi as mpi + +The last two lines load the modules for the construction of the +:ref:`CTHYB solver `. + +Initializing SumkDFT +-------------------- + +We define some parameters, which should be self-explanatory:: + + dft_filename = 'SrVO3' # filename + U = 4.0 # interaction parameters + J = 0.65 + beta = 40 # inverse temperature + loops = 15 # number of DMFT loops + mix = 0.8 # mixing factor of Sigma after solution of the AIM + dc_type = 1 # DC type: 0 FLL, 1 Held, 2 AMF + use_blocks = True # use bloc structure from DFT input + prec_mu = 0.0001 # precision of chemical potential + + +And next, we can initialize the :class:`SumkDFT ` class:: + + SK = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks) + +Initializing the solver +----------------------- + +We also have to specify the :ref:`CTHYB solver ` related settings. +We assume that the DMFT script for SrVO3 is executed on 16 cores. A sufficient set +of parameters for a first guess is:: + + p = {} + # solver + p["random_seed"] = 123 * mpi.rank + 567 + p["length_cycle"] = 200 + p["n_warmup_cycles"] = 100000 + p["n_cycles"] = 1000000 + # tail fit + p["perform_tail_fit"] = True + p["fit_max_moment"] = 4 + p["fit_min_n"] = 30 + p["fit_max_n"] = 60 + +Here we use a tail fit to deal with numerical noise of higher Matsubara frequencies. +For other options and more details on the solver parameters, we refer the user to +the :ref:`CTHYB solver ` documentation. +It is important to note that the solver parameters have to be adjusted for +each material individually. A guide on how to set the tail fit parameters is given +:ref:`below `. + + +The next step is to initialize the +:class:`solver class `. +It consist of two parts: + +#. Calculating the multi-band interaction matrix, and constructing the + interaction Hamiltonian. +#. Initializing the solver class itself. + +The first step is done using methods of the :ref:`TRIQS ` library:: + + n_orb = SK.corr_shells[0]['dim'] + l = SK.corr_shells[0]['l'] + spin_names = ["up","down"] + orb_names = [i for i in range(n_orb)] + # Use GF structure determined by DFT blocks: + gf_struct = [(block, indices) for block, indices in SK.gf_struct_solver[0].iteritems()] + # Construct U matrix for density-density calculations: + Umat, Upmat = U_matrix_kanamori(n_orb=n_orb, U_int=U, J_hund=J) + +We assumed here that we want to use an interaction matrix with +Kanamori definitions of :math:`U` and :math:`J`. + +Next, we construct the Hamiltonian and the solver:: + + h_int = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat) + S = Solver(beta=beta, gf_struct=gf_struct) + +As you see, we take only density-density interactions into +account. Other Hamiltonians with, e.g. with full rotational invariant interactions are: + +* h_int_kanamori +* h_int_slater + +For other choices of the interaction matrices (e.g Slater representation) or +Hamiltonians, we refer to the reference manual of the :ref:`TRIQS ` +library. + +DMFT cycle +---------- + +Now we can go to the definition of the self-consistency step. It consists again +of the basic steps discussed in the :ref:`previous section `, with +some additional refinements:: + + for iteration_number in range(1,loops+1): + if mpi.is_master_node(): print "Iteration = ", iteration_number + + SK.symm_deg_gf(S.Sigma_iw,orb=0) # symmetrizing Sigma + SK.set_Sigma([ S.Sigma_iw ]) # put Sigma into the SumK class + chemical_potential = SK.calc_mu( precision = prec_mu ) # find the chemical potential for given density + S.G_iw << SK.extract_G_loc()[0] # calc the local Green function + mpi.report("Total charge of Gloc : %.6f"%S.G_iw.total_density()) + + # Init the DC term and the real part of Sigma, if no previous runs found: + if (iteration_number==1 and previous_present==False): + dm = S.G_iw.density() + SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) + S.Sigma_iw << SK.dc_imp[0]['up'][0,0] + + # Calculate new G0_iw to input into the solver: + S.G0_iw << S.Sigma_iw + inverse(S.G_iw) + S.G0_iw << inverse(S.G0_iw) + + # Solve the impurity problem: + S.solve(h_int=h_int, **p) + + # Solved. Now do post-solution stuff: + mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density()) + + # Now mix Sigma and G with factor mix, if wanted: + if (iteration_number>1 or previous_present): + if mpi.is_master_node(): + with HDFArchive(dft_filename+'.h5','r') as ar: + mpi.report("Mixing Sigma and G with factor %s"%mix) + S.Sigma_iw << mix * S.Sigma_iw + (1.0-mix) * ar['dmft_output']['Sigma_iw'] + S.G_iw << mix * S.G_iw + (1.0-mix) * ar['dmft_output']['G_iw'] + S.G_iw << mpi.bcast(S.G_iw) + S.Sigma_iw << mpi.bcast(S.Sigma_iw) + + # Write the final Sigma and G to the hdf5 archive: + if mpi.is_master_node(): + with HDFArchive(dft_filename+'.h5','a') as ar: + ar['dmft_output']['iterations'] = iteration_number + ar['dmft_output']['G_0'] = S.G0_iw + ar['dmft_output']['G_tau'] = S.G_tau + ar['dmft_output']['G_iw'] = S.G_iw + ar['dmft_output']['Sigma_iw'] = S.Sigma_iw + + # Set the new double counting: + dm = S.G_iw.density() # compute the density matrix of the impurity problem + SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) + + # Save stuff into the user_data group of hdf5 archive in case of rerun: + SK.save(['chemical_potential','dc_imp','dc_energ']) + + +This is all we need for the DFT+DMFT calculation. +You can see in this code snippet, that all results of this calculation +will be stored in a separate subgroup in the hdf5 file, called `dmft_output`. +Note that this script performs 15 DMFT cycles, but does not check for +convergence. Of course, it would be possible to build in convergence criteria. +A simple check for convergence can be also done if you store multiple quantities +of each iteration and analyse the convergence by hand. In general, it is advisable +to start with a lower statistics (less measurements), but then increase it at a +point close to converged results (e.g. after a few initial iterations). This helps +to keep computational costs low during the first iterations. + +Using the Kanamori Hamiltonian and the parameters above (but on 16 cores), +your self energy after the **first iteration** should look like the +self energy shown below. + +.. image:: images_scripts/SrVO3_Sigma_iw_it1.png + :width: 700 + :align: center + +.. _tailfit: + +Tail fit parameters +------------------- + +A good way to identify suitable tail fit parameters is by "human inspection". +Therefore disabled the tail fitting first:: + + p["perform_tail_fit"] = False + +and perform only one DMFT iteration. The resulting self energy can be tail fitted by hand:: + + Sigma_iw_fit = S.Sigma_iw.copy() + Sigma_iw_fit << tail_fit(S.Sigma_iw, fit_max_moment = 4, fit_min_n = 40, fit_max_n = 160)[0] + +Plot the self energy and adjust the tail fit parameters such that you obtain a +proper fit. The :meth:`fit_tail function ` is part +of the :ref:`TRIQS ` library. + +For a self energy which is going to zero for :math:`i\omega \rightarrow 0` our suggestion is +to start the tail fit (:emphasis:`fit_min_n`) at a Matsubara frequency considerable above the minimum +of the self energy and to stop (:emphasis:`fit_max_n`) before the noise fully takes over. +If it is difficult to find a reasonable fit in this region you should increase +your statistics (number of measurements). Keep in mind that :emphasis:`fit_min_n` +and :emphasis:`fit_max_n` also depend on :math:`\beta`.