nonbonded.py 9.27 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
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
from shutil import copyfile
import os, sys
import numpy as np

class NonbondedScheme():
    """ Abstract class for writing nonbonded interactions """

    def __init__(self, typesA=None, typesB=None, resolution=0.1, rMin=0):
        """If typesA is None, and typesB is None, then """
        self.resolution = resolution
        self.rMin = rMin

    def add_sim_system(self, simSystem):
        self.rMax = simSystem.cutoff
        self.r = np.arange(rMin,rMax,resolution)

    def potential(self, r, typeA, typeB):
        raise NotImplementedError
    
    def write_file(self, filename, typeA, typeB, rMax):
        r = np.arange(self.rMin, rMax+self.resolution, self.resolution)
        u = self.potential(r, typeA, typeB)
        np.savetxt(filename, np.array([r,u]).T)


class LennardJones(NonbondedScheme):
    def potential(self, r, typeA, typeB):
        epsilon = sqrt( typeA.epsilon**2 + typeB.epsilon**2 )
        r0 = 0.5 * (typeA.radius + typeB.radius)
        r6 = (r0/r)**6
        r12 = r6**2
        u = epsilon * (r12-2*r6)
        u[0] = u[1]             # Remove NaN
        return u
LennardJones = LennardJones()

class HalfHarmonic(NonbondedScheme):
    def potential(self, r, typeA, typeB):
        k = 10                   # kcal/mol AA**2
        r0 = (typeA.radius + typeB.radius)
        u =  0.5 * k * (r-r0)**2
        u[r > r0] = np.zeros( np.shape(u[r > r0]) )
        return u
HalfHarmonic = HalfHarmonic()

class TabulatedPotential(NonbondedScheme):
    def __init__(self, tableFile, typesA=None, typesB=None, resolution=0.1, rMin=0):
        """If typesA is None, and typesB is None, then """
        self.tableFile = tableFile
        # self.resolution = resolution
        # self.rMin = rMin

        ## TODO: check that tableFile exists and is regular file
        
    def write_file(self, filename, typeA, typeB, rMax):
        if filename != self.tableFile:
            copyfile(self.tableFile, filename)

## Bonded potentials
60
61
class BasePotential():
    def __init__(self, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
62
63
        self.r0 = r0
        self.rRange = rRange
64
        self.resolution = resolution
65
66
67
68
        self.maxForce = maxForce
        self.prefix = prefix
        self.periodic = False
        self.type_ = "None"
cmaffeo2's avatar
cmaffeo2 committed
69
        self.max_potential = max_potential
70
        self.kscale_ = None     # only used for 
71
72

    def filename(self):
73
74
75
76
        raise NotImplementedError("Not implemented")

    def potential(self,dr):
        raise NotImplementedError("Not implemented")
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91

    def __str__(self):
        return self.filename()

    def write_file(self):
        r = np.arange( self.rRange[0], 
                       self.rRange[1]+self.resolution, 
                       self.resolution )
        dr = r-self.r0

        if self.periodic == True:
            rSpan = self.rRange[1]-self.rRange[0]
            assert(rSpan > 0)
            dr = np.mod( dr+0.5*rSpan, rSpan) - 0.5*rSpan 

92
        u = self.potential(dr)
93

94
95
96
        maxForce = self.maxForce
        if maxForce is not None:
            assert(maxForce > 0)
97
98
99
100
101
            f = np.diff(u)/np.diff(r)
            f[f>maxForce] = maxForce
            f[f<-maxForce] = -maxForce
            u[0] = 0
            u[1:] = np.cumsum(f*np.diff(r))
cmaffeo2's avatar
cmaffeo2 committed
102
103
104
105
106
107
108
109
110
111
112
113
114
115

        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))
            u[0] = 0
            u[1:] = np.cumsum(f*np.diff(r))

        u = u - np.min(u)

116
117
        np.savetxt( self.filename(), np.array([r, u]).T, fmt="%f" )

118
class HarmonicPotential(BasePotential):
119
    def __init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix, correct_geometry=False):
120
121
122
123
124
125
126
127
128
129
        self.k = k
        self.kscale_ = None
        BasePotential.__init__(self, r0, rRange, resolution, maxForce, max_potential, prefix)

    def filename(self):
        return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
                                       self.k*self.kscale_, self.r0)
    def potential(self,dr):
        return 0.5*self.k*dr**2

130
131
    def __hash__(self):
        assert(self.type_ != "None")
cmaffeo2's avatar
cmaffeo2 committed
132
        return hash((self.type_, self.k, self.r0, self.rRange, self.resolution, self.maxForce, self.max_potential, self.prefix, self.periodic))
133
134

    def __eq__(self, other):
cmaffeo2's avatar
cmaffeo2 committed
135
        for a in ("type_", "k", "r0", "rRange", "resolution", "maxForce", "max_potential", "prefix", "periodic"):
136
137
138
139
140
141
142
143
144
145
            if self.__dict__[a] != other.__dict__[a]:
                return False
        return True

# class NonBonded(HarmonicPotential):
#     def _init_hook(self):
#         self.type = "nonbonded"
#         self.kscale_ = 1.0

class HarmonicBond(HarmonicPotential):
146
    def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/", correct_geometry=False, temperature=295):
cmaffeo2's avatar
cmaffeo2 committed
147
        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix)
148
        self.type_ = "gbond" if correct_geometry else "bond"
149
        self.kscale_ = 1.0
150
151
152
153
154
155
        self.correct_geometry = correct_geometry
        self.temperature = temperature

    def potential(self,dr):
        u = HarmonicPotential.potential(self,dr)
        if self.correct_geometry:
156
157
            with np.errstate(divide='ignore',invalid='ignore'):
                du = 2*0.58622592*np.log(dr+self.r0) * self.temperature/295
158
159
160
161
            du[np.logical_not(np.isfinite(du))] = 0
            u = u-du
        return u

162
163

class HarmonicAngle(HarmonicPotential):
cmaffeo2's avatar
cmaffeo2 committed
164
165
    def __init__(self, k, r0, rRange=(0,181), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix)
166
        self.type_ = "angle"
167
        self.kscale_ = (180.0/np.pi)**2
168
169

class HarmonicDihedral(HarmonicPotential):
cmaffeo2's avatar
cmaffeo2 committed
170
171
    def __init__(self, k, r0, rRange=(-180,180), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix)
172
173
        self.periodic = True
        self.type_ = "dihedral"
174
        self.kscale_ = (180.0/np.pi)**2
175

176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
class WLCSKBond(BasePotential):
    """ ## https://aip.scitation.org/doi/full/10.1063/1.4968020 """
    def __init__(self, d, lp, kT, rRange=(0,50), resolution=0.02, maxForce=100, max_potential=None, prefix="potentials/"):
        BasePotential.__init__(self, 0, rRange, resolution, maxForce, max_potential, prefix)
        self.type_ = "wlcbond"
        self.d = d          # separation
        self.lp = lp            # persistence length
        self.kT = kT

    def filename(self):
        return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
                                       self.d, self.lp)
    def potential(self, dr):
        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) )
        max_force = np.diff(u[q2<1][-2:]) / np.diff(dr).mean()
        max_u = u[q2<1][-1]
        max_dr = dr[q2<1][-2]
        assert( max_force >= 0 )
        u[q2>=1] = (dr[q2>=1]-max_dr)*max_force + max_u
        return u

    def __hash__(self):
        assert(self.type_ != "None")
        return hash((self.type_, self.d, self.lp, self.kT, self.rRange, self.resolution, self.maxForce, self.max_potential, self.prefix, self.periodic))

    def __eq__(self, other):
        for a in ("type_", "d", "lp", "kT" "rRange", "resolution", "maxForce", "max_potential", "prefix", "periodic"):
            if self.__dict__[a] != other.__dict__[a]:
                return False
        return True

class WLCSKAngle(BasePotential):
    ## https://aip.scitation.org/doi/full/10.1063/1.4968020
    def __init__(self, d, lp, kT, rRange=(0,181), resolution=0.5, maxForce=None, max_potential=None, prefix="potentials/"):
        BasePotential.__init__(self, 180, rRange, resolution, maxForce, max_potential, prefix)
        self.type_ = "wlcangle"
        self.d = d          # separation
        self.lp = lp            # persistence length
        self.kT = kT

    def filename(self):
        return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
                                       self.d, self.lp)
    def potential(self,dr):
        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

    def __hash__(self):
        assert(self.type_ != "None")
        return hash((self.type_, self.d, self.lp, self.kT, self.rRange, self.resolution, self.maxForce, self.max_potential, self.prefix, self.periodic))

    def __eq__(self, other):
        for a in ("type_", "d", "lp", "kT" "rRange", "resolution", "maxForce", "max_potential", "prefix", "periodic"):
            if self.__dict__[a] != other.__dict__[a]:
                return False
        return True