Extension to infinite lag time may not be a good idea

[YongliangOu]

Hello,

I came across the paper [1] and found this code. In the paper, the diffusion coefficient is extended to infinite lag time. I have doubts about it.

If I understand the theory correctly, the likelihood is based on the two matrixes, one is the rate matrix (with unknown $v$ and $w$ [2]) and the other is transition matrix (from MD trajectories). In paper [1], the likelihood is calculated according to Eq. (33), in which the sum runs over all possible transitions $i\leftarrow j$. However, for the rate matrix $R$, only the diagonal and the secondary diagonal elements are defined while all others are zero.

def init_rate_matrix_pbc(n,v,w):
    """initialize rate matrix from potential vector v and diffusion
    vector w = log(D(i)/delta^2)"""
    assert len(v) == n  # number of bins
    assert len(w) == n
    rate = np.float64(np.zeros((n,n)))  # high precision

    # off-diagonal elements
    diffv = v[1:]-v[:-1] #length n-1  # diffv[i] = v[i+1]-v[i]
    exp1 = w[:n-1]-0.5*diffv
    exp2 = w[:n-1]+0.5*diffv
    rate.ravel()[n::n+1] = np.exp(exp1)[:n-1]
    rate.ravel()[1::n+1] = np.exp(exp2)[:n-1]
    #this amounts to doing:
    #for i in range(n-1):
    #    rate[i+1,i] = np.exp(w[i]-0.5*(v[i+1]-v[i]))
    #    rate[i,i+1] = np.exp(w[i]-0.5*(v[i]-v[i+1]))

    # corners    # periodic boundary conditions
    rate[0,-1]  = np.exp(w[-1]-0.5*(v[0]-v[-1]))
    rate[-1,0]  = np.exp(w[-1]-0.5*(v[-1]-v[0]))
    rate[0,0]   = - rate[1,0] - rate[-1,0]
    rate[-1,-1] = - rate[-2,-1] - rate[0,-1]

    # diagonal elements
    for i in range(1,n-1):
        rate[i,i] = - rate[i-1,i] - rate[i+1,i]
    return rate

For the transition matrix $N$, it records the number of transitions between bins $i\leftarrow j$ with the pre-setted lag time. With larger lag time, diffusion can happen multiple times and, thus, not limited to the neighboring bins. In other words, off-diagonal elements are tended to be filled, while the diagonal and the secondary diagonal elements are tended to be the initial value (zero). This can be verified by the given example files (examples/create-transition-matrix/pbctrans).

# read trajectory
#-----------------
z = read_coor(indir+"ucomo2.dat",rv=True,axis=2) # com=True
# axis = 0,1,2 refers to x,y,z coordinate
# rv = whether the ucomo2.dat file format is used, see below.
# com = True (default) means that the COM-of-lipid columns are skipped

# In case you have a different format, you could use your own
# reading function here.
# This example trajectory only has 200 timesteps, very short.

# construct transition matrix for different lag times
#----------------------------------------------------
for shift in [1,2,3,4,5,10,15,20,30,40,50]:
    lt = shift*dt  # lag time in ps
    print("*"*20+"\nshift "+str(shift)+" lagtime "+str(lt)+" ps")

    A = np.zeros((nbins+2,nbins+2),int)
    print("nbins, len(edges), shape(A), z.shape")
    print(nbins, len(edges), A.shape, z.shape)

    # the N centers of mass of the N O2 molecules,
    # or the 2*N atoms of the N O2 molecules
    for i in range(z.shape[1]):
        traj = z[:,i]-zpbc*np.floor(z[:,i]/zpbc+0.5) # put z=0 in the center
        A = transition_matrix_add1(A,traj,edges,shift=shift)

    count = "pbc"  
    Tmatfile = "%s/transitions.nbins%i.%i.%s.dat"%(outdir,nbins,shift,count)
    write_Tmat_square(A[1:-1,1:-1],Tmatfile,lt,count,edges=edges,dt=dt,dn=shift)

def transition_matrix_add1(A,x,edges,shift=1):
    assert len(x.shape) == 1
    assert shift < len(x)
    assert len(A) == len(edges)+1
    digitized = np.digitize(x,edges)
    #if False:  # pbc
    #  for i in digitized:
    #    assert i > 0
    #    assert i < len(edges)
    #print("min,max,A.shape")
    #print(min(digitized),max(digitized),A.shape)
    #print("min,max",min(x),max(x))
    # periodic boundary conditions: just checking
    print(("check boundary", sum(A[0,:]), sum(A[-1,:]), sum(A[:,0]), sum(A[:,-1])))

    for start,end in zip(digitized[:-shift],digitized[shift:]):
        A[end,start]+=1
    return A

Based on Eq. (33) in paper [1]:

It seems that if we have too many off-diagonal elements in the transition matrix, the first neighbor discretization concept will not work. Imagine for a simulation with long lag time, the diagonal, and the secondary diagonal elements of $N_{ij}$ can all be zeros. Due to the fact that only the diagonal, and the secondary diagonal elements in $R$ are non-zero, the calculated likelihood $L(M)$ will anyway be a fixed number. As a result, the fitted $v$ and $w$ have no physical meaning. Based on it, I suggest that the extension to infinite lag time to determine the final $D$ or $F$ may not be a good idea. People need to pay attention to the transition matrix (by adjusting lag time and bin size such that transitions happen between neighboring bins) before using this code.

[1] https://doi.org/10.1021/acs.jctc.7b00039 [2]

# v -- equal to F, in kBT
# w -- log( D/(dx^2/dt) ), has no unit

Best regards, Yongliang Ou

[annekegh]

Dear Yongliang Ou,

The rate matrix (R) and the matrix that counts transitions (N) do not have the same shape. The rate matrix comes from the diagonalization of the Smoluchowski equation. The matrix that counts transitions can be derived from molecular dynamics trajectories. R is tridiagonal (except for corner elements) because of the discretization scheme. N does not need to be tridiagonal.

To read about the difference between R and N, I suggest going through the paper by Gerhard Hummer, New Journal of Physics, 2005.

By the way, another interesting matrix to look at is the propagator, which is the matrix exponential of R*t. The propagator does not need to be diagonal either.

Kind regards, An

[YongliangOu]

Thanks! I got your point.