Different Persistence Length results when using different time windows.

[acabezon-nextmol]

Hello MDAnalysis community!

I was trying to compute the persistence length for a polymer melt formed by chains of poly(lactic-co-glycolic) acid. My system has 9 chains of PLGA each formed by 156 monomer units. I followed the simulation protocol stated in the paper I was following (paper), which ends up in a 20 ns production run. I've double checked energy and density to make sure the system was properly equilibrated.

Using the generated trajectory I've calculated the persistence length using MDAnalysis and I've found that for the same system I get very different results when using different time windows of the simulation as input. I've run the code for the whole trajectory (20 ns), for the last 10 ns, for the last 5 ns, and for the last 1 ns. The persistence length returned by MDAnalysis (in Angstroms) was 1.859, 4.060, 14.073, and 15.571 respectively. I'm new to polymers so I don't know if this is normal, I was expecting similar values for all the cases, or at least, no such a big difference. Since it was strange for me I've adapted my code to return also the average bond length, the bond autocorrelation, the fit, and the length of the decorrelation. Then I've plotted the bond autocorrelation against the length of the decorrelation. I was expecting the bond autocorrelation to start at 1.0, but the bigger the time window the lower the bond autocorrelation of the first point. As I said I have no clue if this is the normal behavior, I'm also lacking some math knowledge here., I'm attaching the plot for the whole trajectory and for the 5 ns one in case that it can be useful.

Below you can find the code. I use a function get_system_backbones which returns the atom indices that form the backbone of the different chains. It is thought to handle different scenarios so it is a bit long to share it, but it basically finds residues that are only attached to one other residue. Then uses NetworkX to compute the shortest path between all the atom pairs in the selected residues and gets the longest of all of them. For this case there are two backbones of the same length for each chain so I'm basically passing them all since I understand that MDAnalysis computes the average of all the passed backbones (am I wrong in this assumption?.
Thanks in advance for the help!

import sys
import argparse
import os
import MDAnalysis as mda
import matplotlib.pyplot as plt
import MDAnalysis.analysis.polymer as pm
import numpy as np
import json

from get_backbones_updated import get_system_backbones
def persistence_length(tpr_file: str, gro_file: str, trr_file: str, verbose: bool =False) -> dict:
    """Function that computes the persistence length of a polymer system

    :param str tpr_file: Topology of the simulation
    :param str gro_file: Final snapshot of the simulation
    :param str trr_file: Trajectory of the simulation
    :param bool verbose: Default False. Set it True to print the result in terminal
    :return: Dictionary with the persistence length and its units
    :rtype: dict
    """

    persistence_length_dict = {
        'Persistence length': 0.0,
        'Units': 'Å'
    }
    results = {
        "Persistence length" : {
            "value" : 0.0,
            "Units" : "Å"
        },
        "Average Bond Length" : {
            "value" : 0.0,
            "Units" : "Å"
        },
        "Bond Autocorrelation" : {
            "value" : None,
            "Units" : None
        },
        "Fit" : {
            "value" : None,
            "Units" : None
        },
        "Length of the decorrelation" : {
            "value" : None,
            "Units" : None
        }
    }
    system = mda.Universe(tpr_file, trr_file)
    polymer = system.atoms
    backbones_selection = get_system_backbones(gro_file, tpr_file)
    # Determine if molecule has more than one backbone of same length
    backbone_array_dimensions = len(np.shape(next(iter(backbones_selection.values()))))
    # If more than one backbone of the same length was found process them all
    if backbone_array_dimensions == 1:
        backbones = [polymer[backbones_selection[x]] for x in backbones_selection.keys()]
    else:
        backbones = []
        for molecule in backbones_selection.keys():
            for backbone in backbones_selection[molecule]:
                backbones.append(polymer[backbone])
    
    try:
        sorted_bb = [pm.sort_backbone(bb) for bb in backbones]
    except ValueError:
        print(f"Error: Selection command '{selection}' led to a non linear "
              f"polymer")
        sys.exit()
    plen = pm.PersistenceLength(sorted_bb)
    plen.run()
    persistence_length_dict['Persistence length'] = plen.results.lp
    results["Persistence length"]["value"] = plen.results.lp
    results["Average Bond Length"]["value"] = plen.results.lb
    results["Bond Autocorrelation"]["value"] = list(plen.results.bond_autocorrelation)
    results["Fit"]["value"] = list(plen.results.fit)
    results["Length of the decorrelation"]["value"] = list(plen.results.x)
    file_name = tpr_file[:-4] # file without extension
    with open(f"{file_name}_unwrapped.json", "w") as f:
        json.dump(results, f)
        f.close()
    if verbose:
        print(persistence_length_dict)
    return persistence_length_dict