DOC: Effective temperature as an example analysis

[hechtprojects]

Problem description:

It would be a nice example to add the following code for the calculation of the effective temperature as an example to the documentation.

# =====================================================================
# Effective Temperature
# =====================================================================
def susceptibility(utraj, ptraj):
    """
    Calculates the susceptibility from an
    unperturbed and perturbed trajectory
    with the same noise realizations as
    described in Refs. [1]_,[2]_,[3]_.
    
    References
    ----------
    .. [1] D. Villamaina, A. Puglisi, and A. Vulpiani,
       The fluctuation-dissipation relation in sub-
       diffusive systems: the case of granular single-
       file diffusion, J. Stat. Mech. Theory Exp. 2008, 
       L 10001 (2008).
       https://doi.org/10.1088/1742-5468/2008/10/L10001
            
    .. [2] É. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco, 
       and F. van Wijand, How Far from Equilibrium is Active
       Matter?, Phys. Rev. Lett. 117, 038103 (2016).
       https://doi.org/10.1103/PhysRevLett.117.038103
            
    .. [3] G. Ciccotti, G. Jacucci, and I. R. McDonald, "Thought-
       experiments" by molecular dynamics, J. Stat. Phys. 21,
       1 (1979). https://doi.org/10.1007/BF01011477
    
    Parameters
    ----------
    utraj : amep.trajectory.ParticleTrajectory
        Unperturbed trajectory.
    ptraj : amep.trajectory.ParticleTrajectory
        Perturbed trajectory.
        
    Returns
    -------
    np.ndarray
        Susceptibility.
    """
    xi = np.zeros(utraj.nframes)
    
    # loop through all frames
    for i in range(utraj.nframes):
        
        # get x coordinate of all particles
        xp = ptraj[i].unwrapped_coords()[:,0]
        xu = utraj[i].unwrapped_coords()[:,0]
        
        # get perturbing force
        fp = ptraj[i].data('f_p')
        
        # calculate susceptibility
        xi[i] = np.mean((xp-xu)/fp)
        
    return xi


def Teff(utrs, ptrs):
    """
    Effective temperature.
    
    Parameters
    ----------
    utrs : list
        List of unperturbed trajectories.
    ptrs : list
        List of perturbed trajectories.
        
    Returns
    -------
    float
        Average effective temperature.
    np.ndarray
        Effective temperature over time.
    np.ndarray
        Times.
    np.ndarray
        Susceptibility averaged over all given
        trajectories.
    np.ndarray
        MSD averaged over all given trajectories.
    
    """
    # average over all dos
    XI = 0
    MSD = 0
    for i in tqdm(range(100)):
        xi  = susceptibility(utrs[i], ptrs[i])
        msd = amep.evaluate.MSD(utrs[i])
        XI += xi
        MSD += msd.frames
    XI = XI/100
    MSD = MSD/100
    
    # get effective temperature
    T = MSD/XI/4
    
    return T[1:].mean(), T, utrs[0].times, XI, MSD

Suggested content improvement:

No response