StaticHessian Module

tdscha.StaticHessian

StaticHessian(ensemble=None, verbose=False, lanczos_input={})

Bases: object

STATIC HESSIAN

This class is for the advanced computation of the static hessian matrix. This exploit the inversion of the auxiliary systems, which allows for including the fourth order contribution exploiting sparse linear algebra to speedup the calculation.

You can either initialize directly the object passing the ensemble with the configurations, or call the init function after the object has been defined.

You can pass all arguments expected for initializing the DynamicalLanczos.Lanczos object as a dictionary called lanczos_input (for example the modes to be excluded, or the running mode).

Check the DynamicalLanczos.Lanczos docs for more details

Source code in tdscha/StaticHessian.py

def __init__(self, ensemble = None, verbose = False, lanczos_input={}):
    """
    STATIC HESSIAN
    ==============

    This class is for the advanced computation of the static hessian matrix.
    This exploit the inversion of the auxiliary systems, which allows for including the
    fourth order contribution exploiting sparse linear algebra to speedup the calculation.

    You can either initialize directly the object passing the ensemble with the configurations,
    or call the init function after the object has been defined.

    You can pass all arguments expected for initializing the DynamicalLanczos.Lanczos object
    as a dictionary called lanczos_input
    (for example the modes to be excluded, or the running mode).

    Check the DynamicalLanczos.Lanczos docs for more details
    """


    # The minimization variables
    self.vector = None
    self.lanczos_input = lanczos_input
    self.lanczos = tdscha.DynamicalLanczos.Lanczos(ensemble, **lanczos_input)
    self.verbose = False
    self.prefix = "hessian_calculation"
    self.preconitioned = True

    if ensemble is not None:
        self.init(ensemble, verbose)

    # Setup the attribute control
    # Every attribute introduce after this line will raise an exception
    self.__total_attributes__ = [item for item in self.__dict__.keys()]
    self.fixed_attributes = True # This must be the last attribute to be setted

load_status(fname)

Load the current vector from a previous calculation

Source code in tdscha/StaticHessian.py

def load_status(self, fname):
    """
    Load the current vector from a previous calculation
    """
    self.lanczos.load_status(fname + ".npz")

    n_modes = self.lanczos.n_modes
    lenv = (self.lanczos.n_modes * (self.lanczos.n_modes + 1)) // 2
    n_g = (n_modes * (n_modes + 1)) // 2
    n_w = (n_modes * (n_modes**2 + 3*n_modes + 2)) // 6
    self.vector = np.zeros( n_g + n_w, dtype = tdscha.DynamicalLanczos.TYPE_DP)

    self.vector[:] = np.loadtxt(fname)

    # Load all the rest from the json file
    with open(fname + ".json", "r") as fp:
        data = json.load(fp)
        for x in data:
            self.__setattr__(x, data[x])

save_status(fname)

Save the current vector to a file.

Source code in tdscha/StaticHessian.py

def save_status(self, fname):
    """
    Save the current vector to a file.
    """
    np.savetxt(fname, self.vector)
    self.lanczos.save_status(fname)

    # Save all the other info in a json file
    selfdict = vars(self)
    new_dict = {}
    IGNORE = ["lanczos", "vector", "fixed_attributes", "__total_attributes__"]

    for x in selfdict:
        if not x in IGNORE:
            new_dict[x] = selfdict[x]


    with open(fname + ".json", "w") as fp:
        json.dump(new_dict, fp)

init(ensemble=None, verbose=True)

Initialize the StaticHessian with a given ensemble

Source code in tdscha/StaticHessian.py

def init(self, ensemble = None, verbose = True):
    """
    Initialize the StaticHessian with a given ensemble

    Parameters
    ----------
        ensemble : sscha.Ensemble.Ensemble
            The object that contains the configurations
        preconditioner : bool
            If true, the initial vector is set up for the preconditioned minimization.
            Be carefull to use it together with a preconditioned algorithm in the run method,
            otherwise you will start from a very far point from convergence.
        verbose : bool
            If true prints the memory occupied for the calculation
    """

    if ensemble is not None:
        self.lanczos = tdscha.DynamicalLanczos.Lanczos(ensemble, **self.lanczos_input)
        self.lanczos.init()

        n_modes = self.lanczos.n_modes
        lenv = (self.lanczos.n_modes * (self.lanczos.n_modes + 1)) // 2
        n_g = (n_modes * (n_modes + 1)) // 2
        n_w = (n_modes * (n_modes**2 + 3*n_modes + 2)) // 6


        self.vector = np.zeros( n_g + n_w, dtype = tdscha.DynamicalLanczos.TYPE_DP)

        # Initialize vector with the initial guess (the SSCHA matrix)
        counter = 0
        for i in range(n_modes):
            if not self.preconitioned:
                self.vector[counter] = 1 / self.lanczos.w[i]**2
            else:
                self.vector[counter] = 1 / self.lanczos.w[i]

            counter += n_modes - i
    else:
        if self.lanczos.N == 0:
            raise ValueError("Error, you must specify an ensemble or initialize the system loading a lanczos algorithm.")

    self.verbose = verbose


    if verbose:
        print("Memory of StaticHessian initialized.")
        # The seven comes from all the auxiliary varialbes necessary in the gradient computation and the CG
        print("     memory requested: {} Gb of RAM per process".format((self.vector.nbytes) * 3 / 1024**3))
        print("     (excluding memory occupied to store the ensemble)")

run(n_steps, save_dir=None, threshold=1e-06, algorithm='cg-prec', extra_options={})

RUN THE HESSIAN MATRIX CALCULATION

This subroutine runs the algorithm that computes the hessian matrix.

After this subroutine finished, the results are stored in the self.Ginv and selfW.W variables. You can retrieve the Hessian matrix as a CC.Phonons.Phonons object with the retrieve_hessian() subroutine.

Source code in tdscha/StaticHessian.py

def run(self, n_steps, save_dir = None, threshold = 1e-6, algorithm = "cg-prec", extra_options = {}):
    """
    RUN THE HESSIAN MATRIX CALCULATION
    ==================================

    This subroutine runs the algorithm that computes the hessian matrix.

    After this subroutine finished, the results are stored in the
    self.Ginv and selfW.W variables.
    You can retrieve the Hessian matrix as a CC.Phonons.Phonons object
    with the retrieve_hessian() subroutine.


    Parameters
    ----------
        n_steps : int
            The maximum number of steps after which the calculation is stopped.
        save_dir : string
            Path to the directory on which, after each step, the file are saved.
        threshold : float
            The threshold below which the residual are converged.
        algorithm : string
            The algorithm used for the inversion. 
            The algorithms are implemented in the sscha.Tools method, look at them for a more specific guide.
                - cg : Conjugate Gradient
                    This is the 'minimal residual method'. Look at sscha.Tools.minimal_residual_method
                - fom : Full orthogonalization Method
                    This is a method based on krilov subspace as Lanczos. It requires the dimension of the Krilov subspace
                    to be specified in extra_options (Krylov_dimension)
        extra_options : dict
            The extra arguments to be passed to the minimizer function, they depend on the algorithm.
    """
    # Prepare the saving directory
    if save_dir is not None:
        if not os.path.exists(save_dir):
            os.makedirs(save_dir)

    # Prepare the A matrix
    lenv = len(self.vector)
    n_modes = self.lanczos.n_modes
    A = scipy.sparse.linalg.LinearOperator((lenv, lenv), matvec = self.apply_L)

    def prec_mult(x):
        return self.apply_L(x, preconditioner = True)
    A_precond_half = scipy.sparse.linalg.LinearOperator((lenv, lenv), matvec = prec_mult)

    #A_real = sscha.Tools.get_matrix_from_sparse_linop(A)
    #A_prec = sscha.Tools.get_matrix_from_sparse_linop(A_precond_half)
    #np.savetxt("A.dat", A_real) 
    #np.savetxt("A_precond.dat", A_prec)

    # Define the function to save the results
    def callback(x, iters):
        if save_dir is not None:
            xtilde = x
            np.savetxt(os.path.join(save_dir, "{}_step{:05d}.dat".format(self.prefix, iters)), xtilde)
        sys.stdout.flush()


    x0 = self.vector.copy()
    b = np.zeros( lenv, dtype = tdscha.DynamicalLanczos.TYPE_DP)

    # Initialize vector of the known terms
    counter = 0
    for i in range(n_modes):
        b[counter] = 1 
        counter += n_modes - i


    t1 = time.time()
    if algorithm.lower() == "cg":
        res = Tools.minimum_residual_algorithm(A, b, x0, max_iters = n_steps, conv_thr = threshold, callback = callback)
    elif algorithm.lower() == "cg-prec":
        res = Tools.minimum_residual_algorithm_precond(A, b, A_precond_half, x0 = x0, max_iters = n_steps, conv_thr = threshold, callback = callback)
    else:
        raise NotImplementedError("Error, algorithm '{}' not implemented.".format(algorithm))
    t2 = time.time()

    self.vector[:] = res

    if save_dir is not None:
        self.save_status( os.path.join(save_dir, self.prefix))


    if self.verbose:
        print()
        print(r"/======================================\\")
        print(r"|                                      |")
        print(r"|    HESSIAN CALCULATION CONVERGED     |")
        print(r"|                                      |")
        print(r"\\======================================/")
        print()
        print()
        totaltime = t2 -t1 
        hours = totaltime // 3600
        totaltime -= hours * 3600
        mins = totaltime // 60
        totaltime -= mins * 60
        print("Total time to converge: {} h {} m {:.2f} s".format(hours, mins, totaltime))
        print()

run_no_mode_mixing(nsteps, save_dir=None, restart_from_file=False)

RUN THE HESSIAN CALCULATION NEGLECTING MODE MIXING

This algorithm computes the Hessian matrix by applying the static approximation on the diagonal elements of the dynamical Green function for each phonon mode of the auxiliary force constant matrix.

This neglects mode mixing introduced by the "bubble" and higher-order phonon-phonon scattering diagrams, but contains anharmonicity non-perturbatively.

TODO: exploit mode symmetries to avoid to repeat calculations for degenerate modes.

Results

hessian : CC.Phonons.Phonons()
    The free energy Hessian, without the mode-mixing approximation.

Source code in tdscha/StaticHessian.py

def run_no_mode_mixing(self, nsteps, save_dir = None, restart_from_file = False):
    """
    RUN THE HESSIAN CALCULATION NEGLECTING MODE MIXING
    ==================================================

    This algorithm computes the Hessian matrix by applying the static approximation on the diagonal
    elements of the dynamical Green function for each phonon mode of the auxiliary force constant matrix.

    This neglects mode mixing introduced by the "bubble" and higher-order phonon-phonon scattering diagrams,
    but contains anharmonicity non-perturbatively.

    TODO: exploit mode symmetries to avoid to repeat calculations for degenerate modes.

    Parameters
    ----------
        nsteps : int
            The number of steps for each Lanczos iterations to converge
        save_dir : string
            Path to the file on which the Lanczos data will be saved for restarting.
        restart_from_file: bool
            If True, we load the data from save_dir and continue from a previous calculation.
            NOT YET IMPLEMENTED

    Results
    -------
        hessian : CC.Phonons.Phonons()
            The free energy Hessian, without the mode-mixing approximation.
    """
    nmodes = self.lanczos.n_modes

    Gw = np.zeros( (nmodes, nmodes), dtype = np.double)


    if not save_dir == None:
        if not os.path.exists(save_dir):
            os.makedirs(save_dir)

    for i in range(nmodes):
        run_lanczos = True
        # Avoid redundancies in simulating degenerate modes.
        if i > 0:
            if np.abs(self.lanczos.w[i-1] - self.lanczos.w[i]) < 1e-10:
                if self.verbose:
                    print("SKIPPING MODE {} (Degeneracy)".format(i+1))
                run_lanczos = False 

        if run_lanczos:
            self.lanczos.prepare_mode(i)

            if self.verbose:
                print()
                string = "RUNNING MODE {} out of {}".format(i + 1, nmodes)
                print(string)
                print("-"*len(string))
                print()

            new_nsteps = nsteps

            if save_dir is not None:
                final_lanc_filename = os.path.join(save_dir, self.prefix + "_M{:d}_conv".format(i + 1))

                if  restart_from_file:
                    # Check if the lanczos needs to be loaded
                    if os.path.exists( final_lanc_filename ):
                        if self.verbose:
                            print("Loading from {}...".format(final_lanc_filename))
                        self.lanczos.load_status(final_lanc_filename)
                    else:
                        all_files = [x for x in os.listdir(save_dir) if x.startswith(self.prefix + "_M{:d}".format(i+1)) and x.endswith(".npz") and "STEP" in x]
                        if len(all_files) > 0:
                            count = [int(x.split(".")[-2].split("STEP")[-1]) for x in all_files]
                            index = np.argmax(count)
                            filename = os.path.join(save_dir, all_files[index])

                            if self.verbose:
                                print("Loading from {}...".format(filename))
                            self.lanczos.load_status(filename)

                    # Get the number of steps
                    if len(self.lanczos.a_coeffs):
                        new_nsteps = nsteps - len(self.lanczos.a_coeffs)


            if new_nsteps > 0:
                self.lanczos.run_FT(new_nsteps, save_dir = save_dir, verbose = self.verbose, prefix = self.prefix + "_M{:d}".format(i+1))

            # Save the final status
            if save_dir is not None:
                self.lanczos.save_status(final_lanc_filename)

            # END if run_lanczos

        # Get the static limit from the dynamical response funciton
        gf0 = self.lanczos.get_green_function_continued_fraction(np.array([0]), use_terminator = False, smearing = 0)[0]
        Gw[i, i] = gf0

    # Retrive the hessian
    W = np.zeros((nmodes, nmodes, nmodes), dtype = np.double)

    # TODO: Apply the correct preconditioning to the vector
    #       this would allow to use the no_mode_mixing as an ansatz for the full algorithm.
    self.preconitioned = False
    self.vector = self.get_vector(Gw, W)

    return self.retrieve_hessian()

retrieve_hessian(noq=False)

Return the Hessian matrix as a CC.Phonons.Phonons object.

Note that you need to run the Hessian calculation (run method), otherwise this method returns the SSCHA dynamical matrix.

If noq is true, it returns a numpy matrix in the supercell. This enable to get the Hessian even if the dynamical matrix of the Lanczos is not initialized, as it occurs when loading from a file.

Source code in tdscha/StaticHessian.py

def retrieve_hessian(self, noq = False):
    """
    Return the Hessian matrix as a CC.Phonons.Phonons object.

    Note that you need to run the Hessian calculation (run method), otherwise this
    method returns the SSCHA dynamical matrix.

    If noq is true, it returns a numpy matrix in the supercell.
    This enable to get the Hessian even if the dynamical matrix of the Lanczos is not initialized,
    as it occurs when loading from a file.
    """

    new_vector = self.vector.copy()

    if self.preconitioned:
       new_vector = self.apply_L(self.vector, preconditioner = True) 

    Ginv = self.get_G_W(new_vector, ignore_W = True)

    G = np.linalg.inv(Ginv)
    dyn = self.lanczos.pols.dot(G.dot(self.lanczos.pols.T))

    dyn[:,:] *= np.outer(np.sqrt(self.lanczos.m), np.sqrt(self.lanczos.m))

    if noq:
        return dyn
    #CC.symmetries.CustomASR(dyn)

    # We copy the q points from the SSCHA dyn
    nq_tot = len(self.lanczos.dyn.q_tot)
    q_points = np.array(self.lanczos.dyn.q_tot)
    uc_structure = self.lanczos.dyn.structure.copy()
    ss_structure = self.lanczos.dyn.structure.generate_supercell(self.lanczos.dyn.GetSupercell())

    # Compute the dynamical matrix
    dynq = CC.Phonons.GetDynQFromFCSupercell(dyn, q_points, uc_structure, ss_structure)

    # Create the CellConstructor Object.
    hessian_matrix = CC.Phonons.Phonons(self.lanczos.dyn.structure, nqirr = nq_tot) #self.lanczos.dyn.Copy()
    for i in range(nq_tot):
        hessian_matrix.dynmats[i] = dynq[i, :, :]
        hessian_matrix.q_tot[i] = q_points[i, :]

    # Adjust the star according to the symmetries
    hessian_matrix.AdjustQStar()

    #CC.symmetries.CustomASR(hessian_matrix.dynmats[0])

    return hessian_matrix

apply_L(vect, preconditioner=False, inverse_preconditioner=False)

Apply the system matrix to the full array.

Source code in tdscha/StaticHessian.py

def apply_L(self, vect, preconditioner = False, inverse_preconditioner = False):
    """
    Apply the system matrix to the full array.
    """
    if self.verbose:
        t1 = time.time()

    if preconditioner and inverse_preconditioner:
        raise ValueError("Error, only one between preconditioner and inverse_preconditioner can be True")

    lenv = self.lanczos.n_modes
    lenv += (self.lanczos.n_modes * (self.lanczos.n_modes + 1)) // 2

    Ginv, W = self.get_G_W(vect)

    if self.verbose:
        t2 = time.time()
        print("Time to transform the vector in the G and W matrix: {} s".format(t2  - t1))

    ti = time.time()
    for i in range(self.lanczos.n_modes):

        vector = np.zeros(lenv, dtype = tdscha.DynamicalLanczos.TYPE_DP)
        vector[:self.lanczos.n_modes] = Ginv[i, :]
        vector[self.lanczos.n_modes:] = _from_matrix_to_symmetric(W[i, :, :])

        # Here the L application (TODO: Here eventual preconditioning)
        self.lanczos.psi = vector

        if not preconditioner and not inverse_preconditioner:
            outv = self.lanczos.apply_L1_static(vector)
            outv += self.lanczos.apply_anharmonic_static()
        else:
            if preconditioner:
                power = 0.5
            if inverse_preconditioner:
                power = -0.5
            outv = self.lanczos.apply_L1_static(vector, inverse = True, power = power)

        Ginv[i, :] = outv[:self.lanczos.n_modes]
        W[i, :, :] = _from_symmetric_to_matrix(outv[self.lanczos.n_modes:], self.lanczos.n_modes)

        if self.verbose:
            tnew = time.time()
            deltat = tnew - ti
            etatime = deltat * (self.lanczos.n_modes -i -1)
            print("Vector {} / {} done (time = {:.2} s | ETA = {:.2f} s)".format(i +1, self.lanczos.n_modes, deltat, etatime))
            ti = tnew

    if self.verbose:
        t3 = time.time()
        print("Total time to apply the L matrix: {} s".format(t3- t2))


    # Reget the final vector
    vect = self.get_vector(Ginv, W)

    if self.verbose:
        t4 = time.time()
        print("Total time to convert to 1D vector: {} s".format(t4-t3))
        sys.stdout.flush()

    return vect

get_G_W(vector, ignore_W=False)

From a vector of the status, return the G and W tensors. If ignore_W is True, only G is returned

Source code in tdscha/StaticHessian.py

def get_G_W(self, vector, ignore_W = False):
    """
    From a vector of the status, return the G and W tensors.
    If ignore_W is True, only G is returned
    """

    n_modes = self.lanczos.n_modes
    # The first part of the vector describes the G
    G = np.zeros((n_modes, n_modes), dtype = np.double)

    counter = 0
    for i in range(n_modes):
        for j in range(i, n_modes):
            G[i, j] = vector[counter]
            G[j, i] = vector[counter]
            counter += 1

    if ignore_W:
        return G

    W = np.zeros((n_modes, n_modes,  n_modes), dtype = np.double)

    for i in range(n_modes):
        for j in range(i, n_modes):
            for k in range(j, n_modes):
                W[i, j, k] = vector[counter]
                W[i, k, j] = vector[counter]
                W[j, i, k] = vector[counter]
                W[j, k, i] = vector[counter]
                W[k, i, j] = vector[counter]
                W[k, j, i] = vector[counter]
                counter += 1

    return G, W

get_vector(G, W)

Get the symmetric vector from G and W. This subroutine takes only the upper diagonal from the tensors.

Source code in tdscha/StaticHessian.py

def get_vector(self, G, W):
    """
    Get the symmetric vector from G and W.
    This subroutine takes only the upper diagonal from the tensors.
    """
    n_modes = self.lanczos.n_modes

    n_g = (n_modes * (n_modes + 1)) // 2
    n_w = (n_modes * (n_modes**2 + 3*n_modes + 2)) // 6

    vector = np.zeros(n_g + n_w, dtype = np.double)

    counter = 0
    for i in range(n_modes):
        for j in range(i, n_modes):
            vector[counter] = G[i, j]
            counter += 1

    assert counter == n_g 

    for i in range(n_modes):
        for j in range(i, n_modes):
            for k in range(j, n_modes):
                vector[counter] = W[i, j, k]
                counter += 1

    assert counter == n_w + n_g

    return vector

StaticHessian Class

tdscha.StaticHessian.StaticHessian(ensemble=None, verbose=False, lanczos_input={})

Bases: object

STATIC HESSIAN

This class is for the advanced computation of the static hessian matrix. This exploit the inversion of the auxiliary systems, which allows for including the fourth order contribution exploiting sparse linear algebra to speedup the calculation.

You can either initialize directly the object passing the ensemble with the configurations, or call the init function after the object has been defined.

You can pass all arguments expected for initializing the DynamicalLanczos.Lanczos object as a dictionary called lanczos_input (for example the modes to be excluded, or the running mode).

Check the DynamicalLanczos.Lanczos docs for more details

Source code in tdscha/StaticHessian.py

def __init__(self, ensemble = None, verbose = False, lanczos_input={}):
    """
    STATIC HESSIAN
    ==============

    This class is for the advanced computation of the static hessian matrix.
    This exploit the inversion of the auxiliary systems, which allows for including the
    fourth order contribution exploiting sparse linear algebra to speedup the calculation.

    You can either initialize directly the object passing the ensemble with the configurations,
    or call the init function after the object has been defined.

    You can pass all arguments expected for initializing the DynamicalLanczos.Lanczos object
    as a dictionary called lanczos_input
    (for example the modes to be excluded, or the running mode).

    Check the DynamicalLanczos.Lanczos docs for more details
    """


    # The minimization variables
    self.vector = None
    self.lanczos_input = lanczos_input
    self.lanczos = tdscha.DynamicalLanczos.Lanczos(ensemble, **lanczos_input)
    self.verbose = False
    self.prefix = "hessian_calculation"
    self.preconitioned = True

    if ensemble is not None:
        self.init(ensemble, verbose)

    # Setup the attribute control
    # Every attribute introduce after this line will raise an exception
    self.__total_attributes__ = [item for item in self.__dict__.keys()]
    self.fixed_attributes = True # This must be the last attribute to be setted

load_status(fname)

Load the current vector from a previous calculation

Source code in tdscha/StaticHessian.py

def load_status(self, fname):
    """
    Load the current vector from a previous calculation
    """
    self.lanczos.load_status(fname + ".npz")

    n_modes = self.lanczos.n_modes
    lenv = (self.lanczos.n_modes * (self.lanczos.n_modes + 1)) // 2
    n_g = (n_modes * (n_modes + 1)) // 2
    n_w = (n_modes * (n_modes**2 + 3*n_modes + 2)) // 6
    self.vector = np.zeros( n_g + n_w, dtype = tdscha.DynamicalLanczos.TYPE_DP)

    self.vector[:] = np.loadtxt(fname)

    # Load all the rest from the json file
    with open(fname + ".json", "r") as fp:
        data = json.load(fp)
        for x in data:
            self.__setattr__(x, data[x])

save_status(fname)

Save the current vector to a file.

Source code in tdscha/StaticHessian.py

def save_status(self, fname):
    """
    Save the current vector to a file.
    """
    np.savetxt(fname, self.vector)
    self.lanczos.save_status(fname)

    # Save all the other info in a json file
    selfdict = vars(self)
    new_dict = {}
    IGNORE = ["lanczos", "vector", "fixed_attributes", "__total_attributes__"]

    for x in selfdict:
        if not x in IGNORE:
            new_dict[x] = selfdict[x]


    with open(fname + ".json", "w") as fp:
        json.dump(new_dict, fp)

init(ensemble=None, verbose=True)

Initialize the StaticHessian with a given ensemble

Source code in tdscha/StaticHessian.py

def init(self, ensemble = None, verbose = True):
    """
    Initialize the StaticHessian with a given ensemble

    Parameters
    ----------
        ensemble : sscha.Ensemble.Ensemble
            The object that contains the configurations
        preconditioner : bool
            If true, the initial vector is set up for the preconditioned minimization.
            Be carefull to use it together with a preconditioned algorithm in the run method,
            otherwise you will start from a very far point from convergence.
        verbose : bool
            If true prints the memory occupied for the calculation
    """

    if ensemble is not None:
        self.lanczos = tdscha.DynamicalLanczos.Lanczos(ensemble, **self.lanczos_input)
        self.lanczos.init()

        n_modes = self.lanczos.n_modes
        lenv = (self.lanczos.n_modes * (self.lanczos.n_modes + 1)) // 2
        n_g = (n_modes * (n_modes + 1)) // 2
        n_w = (n_modes * (n_modes**2 + 3*n_modes + 2)) // 6


        self.vector = np.zeros( n_g + n_w, dtype = tdscha.DynamicalLanczos.TYPE_DP)

        # Initialize vector with the initial guess (the SSCHA matrix)
        counter = 0
        for i in range(n_modes):
            if not self.preconitioned:
                self.vector[counter] = 1 / self.lanczos.w[i]**2
            else:
                self.vector[counter] = 1 / self.lanczos.w[i]

            counter += n_modes - i
    else:
        if self.lanczos.N == 0:
            raise ValueError("Error, you must specify an ensemble or initialize the system loading a lanczos algorithm.")

    self.verbose = verbose


    if verbose:
        print("Memory of StaticHessian initialized.")
        # The seven comes from all the auxiliary varialbes necessary in the gradient computation and the CG
        print("     memory requested: {} Gb of RAM per process".format((self.vector.nbytes) * 3 / 1024**3))
        print("     (excluding memory occupied to store the ensemble)")

run(n_steps, save_dir=None, threshold=1e-06, algorithm='cg-prec', extra_options={})

RUN THE HESSIAN MATRIX CALCULATION

This subroutine runs the algorithm that computes the hessian matrix.

After this subroutine finished, the results are stored in the self.Ginv and selfW.W variables. You can retrieve the Hessian matrix as a CC.Phonons.Phonons object with the retrieve_hessian() subroutine.

Source code in tdscha/StaticHessian.py

def run(self, n_steps, save_dir = None, threshold = 1e-6, algorithm = "cg-prec", extra_options = {}):
    """
    RUN THE HESSIAN MATRIX CALCULATION
    ==================================

    This subroutine runs the algorithm that computes the hessian matrix.

    After this subroutine finished, the results are stored in the
    self.Ginv and selfW.W variables.
    You can retrieve the Hessian matrix as a CC.Phonons.Phonons object
    with the retrieve_hessian() subroutine.


    Parameters
    ----------
        n_steps : int
            The maximum number of steps after which the calculation is stopped.
        save_dir : string
            Path to the directory on which, after each step, the file are saved.
        threshold : float
            The threshold below which the residual are converged.
        algorithm : string
            The algorithm used for the inversion. 
            The algorithms are implemented in the sscha.Tools method, look at them for a more specific guide.
                - cg : Conjugate Gradient
                    This is the 'minimal residual method'. Look at sscha.Tools.minimal_residual_method
                - fom : Full orthogonalization Method
                    This is a method based on krilov subspace as Lanczos. It requires the dimension of the Krilov subspace
                    to be specified in extra_options (Krylov_dimension)
        extra_options : dict
            The extra arguments to be passed to the minimizer function, they depend on the algorithm.
    """
    # Prepare the saving directory
    if save_dir is not None:
        if not os.path.exists(save_dir):
            os.makedirs(save_dir)

    # Prepare the A matrix
    lenv = len(self.vector)
    n_modes = self.lanczos.n_modes
    A = scipy.sparse.linalg.LinearOperator((lenv, lenv), matvec = self.apply_L)

    def prec_mult(x):
        return self.apply_L(x, preconditioner = True)
    A_precond_half = scipy.sparse.linalg.LinearOperator((lenv, lenv), matvec = prec_mult)

    #A_real = sscha.Tools.get_matrix_from_sparse_linop(A)
    #A_prec = sscha.Tools.get_matrix_from_sparse_linop(A_precond_half)
    #np.savetxt("A.dat", A_real) 
    #np.savetxt("A_precond.dat", A_prec)

    # Define the function to save the results
    def callback(x, iters):
        if save_dir is not None:
            xtilde = x
            np.savetxt(os.path.join(save_dir, "{}_step{:05d}.dat".format(self.prefix, iters)), xtilde)
        sys.stdout.flush()


    x0 = self.vector.copy()
    b = np.zeros( lenv, dtype = tdscha.DynamicalLanczos.TYPE_DP)

    # Initialize vector of the known terms
    counter = 0
    for i in range(n_modes):
        b[counter] = 1 
        counter += n_modes - i


    t1 = time.time()
    if algorithm.lower() == "cg":
        res = Tools.minimum_residual_algorithm(A, b, x0, max_iters = n_steps, conv_thr = threshold, callback = callback)
    elif algorithm.lower() == "cg-prec":
        res = Tools.minimum_residual_algorithm_precond(A, b, A_precond_half, x0 = x0, max_iters = n_steps, conv_thr = threshold, callback = callback)
    else:
        raise NotImplementedError("Error, algorithm '{}' not implemented.".format(algorithm))
    t2 = time.time()

    self.vector[:] = res

    if save_dir is not None:
        self.save_status( os.path.join(save_dir, self.prefix))


    if self.verbose:
        print()
        print(r"/======================================\\")
        print(r"|                                      |")
        print(r"|    HESSIAN CALCULATION CONVERGED     |")
        print(r"|                                      |")
        print(r"\\======================================/")
        print()
        print()
        totaltime = t2 -t1 
        hours = totaltime // 3600
        totaltime -= hours * 3600
        mins = totaltime // 60
        totaltime -= mins * 60
        print("Total time to converge: {} h {} m {:.2f} s".format(hours, mins, totaltime))
        print()

run_no_mode_mixing(nsteps, save_dir=None, restart_from_file=False)

RUN THE HESSIAN CALCULATION NEGLECTING MODE MIXING

This algorithm computes the Hessian matrix by applying the static approximation on the diagonal elements of the dynamical Green function for each phonon mode of the auxiliary force constant matrix.

This neglects mode mixing introduced by the "bubble" and higher-order phonon-phonon scattering diagrams, but contains anharmonicity non-perturbatively.

TODO: exploit mode symmetries to avoid to repeat calculations for degenerate modes.

Results

hessian : CC.Phonons.Phonons()
    The free energy Hessian, without the mode-mixing approximation.

Source code in tdscha/StaticHessian.py

def run_no_mode_mixing(self, nsteps, save_dir = None, restart_from_file = False):
    """
    RUN THE HESSIAN CALCULATION NEGLECTING MODE MIXING
    ==================================================

    This algorithm computes the Hessian matrix by applying the static approximation on the diagonal
    elements of the dynamical Green function for each phonon mode of the auxiliary force constant matrix.

    This neglects mode mixing introduced by the "bubble" and higher-order phonon-phonon scattering diagrams,
    but contains anharmonicity non-perturbatively.

    TODO: exploit mode symmetries to avoid to repeat calculations for degenerate modes.

    Parameters
    ----------
        nsteps : int
            The number of steps for each Lanczos iterations to converge
        save_dir : string
            Path to the file on which the Lanczos data will be saved for restarting.
        restart_from_file: bool
            If True, we load the data from save_dir and continue from a previous calculation.
            NOT YET IMPLEMENTED

    Results
    -------
        hessian : CC.Phonons.Phonons()
            The free energy Hessian, without the mode-mixing approximation.
    """
    nmodes = self.lanczos.n_modes

    Gw = np.zeros( (nmodes, nmodes), dtype = np.double)


    if not save_dir == None:
        if not os.path.exists(save_dir):
            os.makedirs(save_dir)

    for i in range(nmodes):
        run_lanczos = True
        # Avoid redundancies in simulating degenerate modes.
        if i > 0:
            if np.abs(self.lanczos.w[i-1] - self.lanczos.w[i]) < 1e-10:
                if self.verbose:
                    print("SKIPPING MODE {} (Degeneracy)".format(i+1))
                run_lanczos = False 

        if run_lanczos:
            self.lanczos.prepare_mode(i)

            if self.verbose:
                print()
                string = "RUNNING MODE {} out of {}".format(i + 1, nmodes)
                print(string)
                print("-"*len(string))
                print()

            new_nsteps = nsteps

            if save_dir is not None:
                final_lanc_filename = os.path.join(save_dir, self.prefix + "_M{:d}_conv".format(i + 1))

                if  restart_from_file:
                    # Check if the lanczos needs to be loaded
                    if os.path.exists( final_lanc_filename ):
                        if self.verbose:
                            print("Loading from {}...".format(final_lanc_filename))
                        self.lanczos.load_status(final_lanc_filename)
                    else:
                        all_files = [x for x in os.listdir(save_dir) if x.startswith(self.prefix + "_M{:d}".format(i+1)) and x.endswith(".npz") and "STEP" in x]
                        if len(all_files) > 0:
                            count = [int(x.split(".")[-2].split("STEP")[-1]) for x in all_files]
                            index = np.argmax(count)
                            filename = os.path.join(save_dir, all_files[index])

                            if self.verbose:
                                print("Loading from {}...".format(filename))
                            self.lanczos.load_status(filename)

                    # Get the number of steps
                    if len(self.lanczos.a_coeffs):
                        new_nsteps = nsteps - len(self.lanczos.a_coeffs)


            if new_nsteps > 0:
                self.lanczos.run_FT(new_nsteps, save_dir = save_dir, verbose = self.verbose, prefix = self.prefix + "_M{:d}".format(i+1))

            # Save the final status
            if save_dir is not None:
                self.lanczos.save_status(final_lanc_filename)

            # END if run_lanczos

        # Get the static limit from the dynamical response funciton
        gf0 = self.lanczos.get_green_function_continued_fraction(np.array([0]), use_terminator = False, smearing = 0)[0]
        Gw[i, i] = gf0

    # Retrive the hessian
    W = np.zeros((nmodes, nmodes, nmodes), dtype = np.double)

    # TODO: Apply the correct preconditioning to the vector
    #       this would allow to use the no_mode_mixing as an ansatz for the full algorithm.
    self.preconitioned = False
    self.vector = self.get_vector(Gw, W)

    return self.retrieve_hessian()

retrieve_hessian(noq=False)

Return the Hessian matrix as a CC.Phonons.Phonons object.

Note that you need to run the Hessian calculation (run method), otherwise this method returns the SSCHA dynamical matrix.

If noq is true, it returns a numpy matrix in the supercell. This enable to get the Hessian even if the dynamical matrix of the Lanczos is not initialized, as it occurs when loading from a file.

Source code in tdscha/StaticHessian.py

def retrieve_hessian(self, noq = False):
    """
    Return the Hessian matrix as a CC.Phonons.Phonons object.

    Note that you need to run the Hessian calculation (run method), otherwise this
    method returns the SSCHA dynamical matrix.

    If noq is true, it returns a numpy matrix in the supercell.
    This enable to get the Hessian even if the dynamical matrix of the Lanczos is not initialized,
    as it occurs when loading from a file.
    """

    new_vector = self.vector.copy()

    if self.preconitioned:
       new_vector = self.apply_L(self.vector, preconditioner = True) 

    Ginv = self.get_G_W(new_vector, ignore_W = True)

    G = np.linalg.inv(Ginv)
    dyn = self.lanczos.pols.dot(G.dot(self.lanczos.pols.T))

    dyn[:,:] *= np.outer(np.sqrt(self.lanczos.m), np.sqrt(self.lanczos.m))

    if noq:
        return dyn
    #CC.symmetries.CustomASR(dyn)

    # We copy the q points from the SSCHA dyn
    nq_tot = len(self.lanczos.dyn.q_tot)
    q_points = np.array(self.lanczos.dyn.q_tot)
    uc_structure = self.lanczos.dyn.structure.copy()
    ss_structure = self.lanczos.dyn.structure.generate_supercell(self.lanczos.dyn.GetSupercell())

    # Compute the dynamical matrix
    dynq = CC.Phonons.GetDynQFromFCSupercell(dyn, q_points, uc_structure, ss_structure)

    # Create the CellConstructor Object.
    hessian_matrix = CC.Phonons.Phonons(self.lanczos.dyn.structure, nqirr = nq_tot) #self.lanczos.dyn.Copy()
    for i in range(nq_tot):
        hessian_matrix.dynmats[i] = dynq[i, :, :]
        hessian_matrix.q_tot[i] = q_points[i, :]

    # Adjust the star according to the symmetries
    hessian_matrix.AdjustQStar()

    #CC.symmetries.CustomASR(hessian_matrix.dynmats[0])

    return hessian_matrix

apply_L(vect, preconditioner=False, inverse_preconditioner=False)

Apply the system matrix to the full array.

Source code in tdscha/StaticHessian.py

def apply_L(self, vect, preconditioner = False, inverse_preconditioner = False):
    """
    Apply the system matrix to the full array.
    """
    if self.verbose:
        t1 = time.time()

    if preconditioner and inverse_preconditioner:
        raise ValueError("Error, only one between preconditioner and inverse_preconditioner can be True")

    lenv = self.lanczos.n_modes
    lenv += (self.lanczos.n_modes * (self.lanczos.n_modes + 1)) // 2

    Ginv, W = self.get_G_W(vect)

    if self.verbose:
        t2 = time.time()
        print("Time to transform the vector in the G and W matrix: {} s".format(t2  - t1))

    ti = time.time()
    for i in range(self.lanczos.n_modes):

        vector = np.zeros(lenv, dtype = tdscha.DynamicalLanczos.TYPE_DP)
        vector[:self.lanczos.n_modes] = Ginv[i, :]
        vector[self.lanczos.n_modes:] = _from_matrix_to_symmetric(W[i, :, :])

        # Here the L application (TODO: Here eventual preconditioning)
        self.lanczos.psi = vector

        if not preconditioner and not inverse_preconditioner:
            outv = self.lanczos.apply_L1_static(vector)
            outv += self.lanczos.apply_anharmonic_static()
        else:
            if preconditioner:
                power = 0.5
            if inverse_preconditioner:
                power = -0.5
            outv = self.lanczos.apply_L1_static(vector, inverse = True, power = power)

        Ginv[i, :] = outv[:self.lanczos.n_modes]
        W[i, :, :] = _from_symmetric_to_matrix(outv[self.lanczos.n_modes:], self.lanczos.n_modes)

        if self.verbose:
            tnew = time.time()
            deltat = tnew - ti
            etatime = deltat * (self.lanczos.n_modes -i -1)
            print("Vector {} / {} done (time = {:.2} s | ETA = {:.2f} s)".format(i +1, self.lanczos.n_modes, deltat, etatime))
            ti = tnew

    if self.verbose:
        t3 = time.time()
        print("Total time to apply the L matrix: {} s".format(t3- t2))


    # Reget the final vector
    vect = self.get_vector(Ginv, W)

    if self.verbose:
        t4 = time.time()
        print("Total time to convert to 1D vector: {} s".format(t4-t3))
        sys.stdout.flush()

    return vect

get_G_W(vector, ignore_W=False)

From a vector of the status, return the G and W tensors. If ignore_W is True, only G is returned

Source code in tdscha/StaticHessian.py

def get_G_W(self, vector, ignore_W = False):
    """
    From a vector of the status, return the G and W tensors.
    If ignore_W is True, only G is returned
    """

    n_modes = self.lanczos.n_modes
    # The first part of the vector describes the G
    G = np.zeros((n_modes, n_modes), dtype = np.double)

    counter = 0
    for i in range(n_modes):
        for j in range(i, n_modes):
            G[i, j] = vector[counter]
            G[j, i] = vector[counter]
            counter += 1

    if ignore_W:
        return G

    W = np.zeros((n_modes, n_modes,  n_modes), dtype = np.double)

    for i in range(n_modes):
        for j in range(i, n_modes):
            for k in range(j, n_modes):
                W[i, j, k] = vector[counter]
                W[i, k, j] = vector[counter]
                W[j, i, k] = vector[counter]
                W[j, k, i] = vector[counter]
                W[k, i, j] = vector[counter]
                W[k, j, i] = vector[counter]
                counter += 1

    return G, W

get_vector(G, W)

Get the symmetric vector from G and W. This subroutine takes only the upper diagonal from the tensors.

Source code in tdscha/StaticHessian.py

def get_vector(self, G, W):
    """
    Get the symmetric vector from G and W.
    This subroutine takes only the upper diagonal from the tensors.
    """
    n_modes = self.lanczos.n_modes

    n_g = (n_modes * (n_modes + 1)) // 2
    n_w = (n_modes * (n_modes**2 + 3*n_modes + 2)) // 6

    vector = np.zeros(n_g + n_w, dtype = np.double)

    counter = 0
    for i in range(n_modes):
        for j in range(i, n_modes):
            vector[counter] = G[i, j]
            counter += 1

    assert counter == n_g 

    for i in range(n_modes):
        for j in range(i, n_modes):
            for k in range(j, n_modes):
                vector[counter] = W[i, j, k]
                counter += 1

    assert counter == n_w + n_g

    return vector