Skip to content

interactions

Module documentation for interactions.

API Reference

AbstractPotential

Abstract class for writing potentials

Source code in mrdna/arbdmodel/submodule/interactions.py
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
class AbstractPotential(metaclass=ABCMeta):
    """ Abstract class for writing potentials """

    def __init__(self, range_=(0,None), resolution=0.1, 
                 max_force = None, max_potential = None, zero='last'):
        self.range_ = range_
        self.resolution = resolution
        self.max_force = max_force
        self.max_potential = max_potential
        self.zero = zero

    @property
    def range_(self):
        return self.__range_
    @range_.setter
    def range_(self,value):
        assert(len(value) == 2)
        self.__range_ = tuple(value)

    @property
    def periodic(self):
        return False

    @abstractmethod
    def potential(self, r, types=None):
        raise NotImplementedError

    def __remove_nans(self, u):
        z = lambda x: np.argwhere(x).flatten()
        nans=np.isnan(u)
        u[nans] = scipy.interpolate.interp1d(z(~nans), u[~nans], fill_value='extrapolate')(z(nans))

    def __remove_nans_bak(self, u):
        left = np.isnan(u[:-1])
        right = np.where(left)[0]+1
        while np.any(np.isnan(u[right])):
            right[np.isnan(u[right])] += 1
        u[:-1][left] = u[right]

    def filename(self, types=None):
        raise NotImplementedError('Inherited potential objects should overload this function')

    def _cap_potential(self, r, u):
        self.__remove_nans(u)

        if self.zero == 'min':
            u = u - np.min(u)
        elif self.zero == 'first':
            u = u - u[0]
        elif self.zero == 'last':
            u = u - u[-1]
        else:
            raise ValueError('Unrecognized option for "zero" argument')

        max_force = self.max_force
        if max_force is not None:
            assert(max_force > 0)
            f = np.diff(u)/np.diff(r)
            f[f>max_force] = max_force
            f[f<-max_force] = -max_force
            u[0] = 0
            u[1:] = np.cumsum(f*np.diff(r))

        if self.max_potential is not None:
            f = np.diff(u)/np.diff(r)
            ids = np.where( 0.5*(u[1:]+u[:-1]) > self.max_potential )[0]

            # w = np.sqrt(2*self.max_potential/self.k)
            # drAvg = 0.5*(np.abs(dr[ids]) + np.abs(dr[ids+1]))
            # f[ids] = f[ids] * np.exp(-(drAvg-w)/(w))
            f[ids] = 0

            u[0] = 0
            u[1:] = np.cumsum(f*np.diff(r))

        if self.zero == 'min':
            u = u - np.min(u)
        elif self.zero == 'first':
            u = u - u[0]
        elif self.zero == 'last':
            u = u - u[-1]
        else:
            raise ValueError('Unrecognized option for "zero" argument')
        return u

    def write_file(self, filename=None, types=None):
        if filename is None:
            filename = self.filename(types)
        rmin,rmax = self.range_
        r = np.arange(rmin, rmax+self.resolution, self.resolution)
        with np.errstate(divide='ignore',invalid='ignore'):
            u = self.potential(r, types)

        u = self._cap_potential(r,u)

        np.savetxt(filename, np.array([r,u]).T)

    def __hash__(self):
        return hash((self.range_, self.resolution, self.max_force, self.max_potential, self.zero))

    def __eq__(self, other):
        for a in ("range_", "resolution", "max_force", "max_potential", "zero"):
            if getattr(self,a) != getattr(other,a):
                return False
        return type(self).__name__ == type(other).__name__

    ## Want the same objects after copy operations
    def __copy__(self):
        return self
    def __deepcopy__(self, memo):
        return self

WLCSKAngle

Bases: WLCSKPotential

https://aip.scitation.org/doi/full/10.1063/1.4968020

Source code in mrdna/arbdmodel/submodule/interactions.py
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
class WLCSKAngle(WLCSKPotential):
    """ ## https://aip.scitation.org/doi/full/10.1063/1.4968020 """
    def __init__(self, d, lp, kT, **kwargs):
        if 'range_' not in kwargs: kwargs['range_'] = (0,181)
        if 'resolution' not in kwargs: kwargs['resolution'] = 0.5
        WLCSKPotential.__init__(self, d, lp, kT, **kwargs)

    @property
    def type_(self):
        return "wlcangle"

    def potential(self, r, types=None):
        dr = r - 180
        nk = self.d / (2*self.lp)
        p1,p2,p3,p4 = -1.237, 0.8105, -1.0243, 0.4595
        C = (1 + p1*(2*nk) + p2*(2*nk)**2) / (2*nk+p3*(2*nk)**2+p4*(2*nk)**3)
        u = self.kT * C * (1-np.cos(dr * np.pi / 180))
        return u

WLCSKBond

Bases: WLCSKPotential

https://aip.scitation.org/doi/full/10.1063/1.4968020

Source code in mrdna/arbdmodel/submodule/interactions.py
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
class WLCSKBond(WLCSKPotential):
    """ ## https://aip.scitation.org/doi/full/10.1063/1.4968020 """
    def __init__(self, d, lp, kT, **kwargs):
        if 'range_' not in kwargs: kwargs['range_'] = (0,50)
        if 'resolution' not in kwargs: kwargs['resolution'] = 0.02
        if 'max_force' not in kwargs: kwargs['max_force'] = 100
        WLCSKPotential.__init__(self, d, lp, kT, **kwargs)

    @property
    def type_(self):
        return "wlcbond"

    def potential(self, r, types=None):
        dr = r
        nk = self.d / (2*self.lp)
        q2 = (dr / self.d)**2
        a1,a2 = 1, -7.0/(2*nk)
        a3 = 3.0/32 - 3.0/(8*nk) - 6.0/(4*nk**2)
        p0,p1,p2,p3,p4 = 13.0/32, 3.4719,2.5064,-1.2906,0.6482
        a4 = (p0 + p1/(2*nk) + p2*(2*nk)**-2) / (1+p3/(2*nk)+p4*(2*nk)**-2)
        with np.errstate(divide='ignore',invalid='ignore'):
            u = self.kT * nk * ( a1/(1-q2) - a2*np.log(1-q2) + a3*q2 - 0.5*a4*q2*(q2-2) )
        return u

WLCSKPotential

Bases: HarmonicBondedPotential

https://aip.scitation.org/doi/full/10.1063/1.4968020

Source code in mrdna/arbdmodel/submodule/interactions.py
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
class WLCSKPotential(HarmonicBondedPotential, metaclass=ABCMeta):
    """ ## https://aip.scitation.org/doi/full/10.1063/1.4968020 """
    def __init__(self, d, lp, kT, **kwargs):
        ## Note, we're leveraging the HarmonicBondedPotential base class and set k to lp here, but it isn't proper
        HarmonicBondedPotential.__init__(self, d, lp, **kwargs)
        self.d = d          # separation
        self.lp = lp            # persistence length
        self.kT = kT

    def filename(self,types=None):
        assert(types is None)
        return "%s%s-%.3f-%.3f.dat" % (self.filename_prefix, self.type_,
                                       self.d, self.lp)

    def __hash__(self):
        return hash((self.d, self.lp, self.kT, HarmonicBondedPotential.__hash__(self)))

    def __eq__(self, other):
        for a in ("d", "lp", "kT"):
            if self.__dict__[a] != other.__dict__[a]:
                return False
        return HarmonicBondedPotential.__eq__(self,other)