Division by zero error

[WaDhondt]

Dear Lukas

I am trying to use PyDHAMed to compute a PMF for a small molecule crossing a membrane bilayer. I performed umbrella sampling in 1.5Å windows, spanning approx. 7nm (57 umbrellas in total).

Following your code, I first created the list of transition matrices for each of the umbrellas:

def create_transition_matrix(traj: npt.NDArray, microstates: npt.NDArray)->npt.NDArray:
    """
    Parameters
    --------------
    traj: np.array of size n
    microstates: np.array of size m

    First, the trajectory is binned into the provided
    microstates. Next, a transition matrix is calculated
    giving the number of transitions between each of the
    microstates.
    
    Returns: m x m matrix with the transitions between 
    all possible pairs of bins in subsequent trajectory frames. 
    """
    traj_binned = np.digitize(traj, microstates)
    n_states = microstates.shape[0]+1 # digitize is 1 indexed.
    mat = np.zeros((n_states, n_states))
    for (x, y), c in six.iteritems(Counter(zip(traj_binned[:-1], traj_binned[1:]))):
        mat[int(y), int(x)] = c
    return mat

colvar_files = glob.glob("../gromacs/plumed_us_*-colvar.dat")

transition_matrices = list()
for i, f in enumerate(colvar_files):
    print("Processing trajectory: ", f)
    traj = np.loadtxt(f)[500:, 1]
    transition_matrices.append(create_transition_matrix(traj, microstates)[1:, 1:])
with open("NFP_membrane_transition_counts.pickle", "wb") as outfile:
    pickle.dump(transition_matrices, outfile)

Next, I calculated the bias grid (in kBT units) as follows:

def create_biasses(q0, kappa, microstates)->npt.NDArray:
    bias = (0.5 * force_constants[:, np.newaxis] * (centers[:, np.newaxis] - microstates)**2).T
    # convert to kBT
    return (bias*1000)/(R*T)

#col0 = bias number
#col1 = bias center position
#col2 = bias force constant

umbrellas = np.loadtxt("../gromacs/wham_input.txt", skiprows=1)
centers = umbrellas[:, 1]
force_constants = umbrellas[:, 2]

biasses = create_biasses(centers, force_constants, microstates)
np.save("NFP_membrane_biasses.npy", biasses)

, where the input file simply contains the centers and force constants of the harmonic potentials. But when I try to use run_dhamedI cannot seem to get sensible results. When I use jit_gradient=False, the calculation finished but I do not get a sensible result (see figure).

When I use jit_gradient=True, I get a division by zero error at '_loop_grad_dhamed_likelihood_0'. Would you have an idea of what is going wrong?

Best,
Warre

[WaDhondt]

I just noticed that I was using the local variables force_constantsand centers in create_biasses instead of the function arguments, but that does not change the outcome.