A python function for create the energy landscape between two conformations using linear interpolation morphing

energyBarrier.py

"""
Mahd Sadjadi (2018)
This is a set of functions to compute the 
energy barrier between two spring networks. 
The morphing is done using the linear interpolation.

Example:
-------
s1 = np.random.rand(10,2)
s2 = np.random.rand(10,2)
b = np.random.randint(low=0,high=10,size=(15,2))
lambdas, energies, morpheds = basinBuilder(s1,s2,b)
"""

import numpy as np

def energy(sites, bonds, L0):
    """
    Compute the energy of a given spring networks
    with free boundary conditions.
    
    Inputs:
    -------
    sites: coordinates of vertices (N*dim)
    bonds: a list of spring ends (M*2)
    L0: rest/natural length of springs (M*1)
    
    Output:
    -------
    energy: energy of the spring network.
    """
    dim = sites1.shape[1]
    dr = np.diff(sites[bonds],axis=1).reshape(-1,dim) # vectors
    L = np.linalg.norm(dr,axis=1)
    energy = 0.5*np.sum((L-L0)**2)
    return energy

def basinBuilder(sites1, sites2, bonds, L0=None, low=-0.5, high=1.5, npoint=100):
    """
    Compute the energy and barrier parameter
    between two given conformations of a spring 
    network.
    with free boundary conditions.
    
    Inputs:
    -------
    sites1: first conformation coordinates (N*dim)
    sites2: second conformation coordinates (N*dim)
    bonds: a list of spring ends (M*2)
    L0: rest/natural length of springs (M*)
        if not provided, the lengths are computed
        from the first conformation.
    low: the lowest value of mo
    hing paramter
    high: the highest value of morphing paramter
    npoint: number of morphings
    
    Output:
    -------
    lambdas: morphing paramter
    morpheds: a list of morphed coordinates
    energy: energy of the spring network.
    """
    if L0 is None:
        dim = sites1.shape[1]
        dr = np.diff(sites1[bonds],axis=1).reshape(-1,dim) # vectors
        L0 = np.linalg.norm(dr,axis=1)
    
    diff = sites2 - sites1
    lambdas = np.linspace(low,high,npoint)
    morpheds = sites1 + np.einsum('i,jk->ijk', lambdas, diff)
    energies = [energy(morphed,bonds,L0) for morphed in morpheds]
    return lambdas, energies, morpheds