Merge pull request #34 from sandialabs/2d

2d

Author: mrbuche

.gitignore

@@ -3,6 +3,6 @@
 **.pyc
 *dist/
 *build/
-*statmechcrack.egg-info/
+**.egg-info/
 **.coverage
 *.vscode/

MiscPlotting/silica.py

@@ -0,0 +1,36 @@
+import numpy as np
+import matplotlib.pyplot as pl
+
+
+wn = 1050*100 #1/m, si-o-si stretching frequency, measured from FTIR
+b = 5e-10   # approximate lattice spacing
+l = 1.55e-10 # Si-O bond length
+ws = 7520 # m/s, elastic wave speed
+E = 87e9 # modulus C11
+udd = 8e-19 #J, Si-O bond energy in water
+
+# set xdd to be a scaling of equilibrium bond length, if xdd=0 there is no Bell term acting
+xdd = 0.0*l
+k = 1.38e-23 #J/K
+T = 293 # K
+
+
+def velSM(K):
+    R = K**2/E
+    f = np.sqrt(R*E*b**3)
+
+    omega = ws*wn
+
+    v = b*omega/np.pi
+    v *= np.exp( (f*xdd-udd)/(k*T) )
+    v *= np.sinh( 0.5*(R*b**2)/(k*T) )
+    return v
+
+npt = 100
+kr = np.linspace(0.1,1.0,npt)
+vr = np.zeros(npt)
+for i in range(npt):
+    vr[i] = velSM(kr[i]*1e6)
+
+pl.semilogy(kr,vr,'k-')
+pl.show()

markovChain.py

@@ -0,0 +1,122 @@
+from multiprocessing import Pool
+import numpy as np
+from statmechcrack import CrackQ2D
+
+def Nn(L,W,Nnom,n=0,verbose=False):
+    if n==0:
+        return  Nnom * np.ones(W, dtype=int) 
+    elif n==1:
+        pass
+    
+
+def gen_N(L,W,Nnom,verbose=False):
+    N0 =  Nnom * np.ones(W, dtype=int) 
+    N = [N0]
+    if verbose:
+        print(N[-1])
+    for i in range(W):
+        N1 = Nnom * np.ones(W, dtype=int) 
+        N1[i] = Nnom+1
+        N.append(N1)
+        if verbose:
+            print(N[-1])
+
+    for i in range(W):
+        N1 = Nnom * np.ones(W, dtype=int) 
+        N1[i] = Nnom+2
+        N.append(N1)
+        if verbose:
+            print(N[-1])
+
+    for i in range(W-1):
+        N1 = Nnom * np.ones(W, dtype=int) 
+        N1[i] = Nnom+1
+        N1[i+1] = Nnom+1
+        N.append(N1)
+        if verbose:
+            print(N[-1])
+
+    for i in range(W-2):
+        N1 = Nnom * np.ones(W, dtype=int) 
+        N1[i] = Nnom+1
+        N1[i+1] = Nnom+1
+        N1[i+2] = Nnom+1
+        N.append(N1)
+        if verbose:
+            print(N[-1])
+
+    for i in range(W-1):
+        N1 = Nnom * np.ones(W, dtype=int) 
+        N1[i] = Nnom+1
+        N1[i+1] = Nnom+2
+        N.append(N1)
+        if verbose:
+            print(N[-1])
+
+    for i in range(W-2):
+        N1 = Nnom * np.ones(W, dtype=int) 
+        N1[i] = Nnom+2
+        N1[i+1] = Nnom+1
+        N1[i+2] = Nnom+1
+        N.append(N1)
+        if verbose:
+            print(N[-1])
+
+    for i in range(W-2):
+        N1 = Nnom * np.ones(W, dtype=int) 
+        N1[i] = Nnom+1
+        N1[i+1] = Nnom+2
+        N1[i+2] = Nnom+1
+        N.append(N1)
+        if verbose:
+            print(N[-1])
+    
+    return N
+
+def calc_rates(L,W,N,verbose=False):
+    model = CrackQ2D(L=L,W=W,N=N)
+    rates = np.zeros(model.W)
+    for k in range(model.W):
+        rates[k] = model.k_isometric(2*np.ones(model.W),k)
+    if verbose:
+        print(N)
+        print(rates)
+    return [N,rates]
+
+def saveRates(N,rates,fname):
+    ofile = open(fname,'w')
+    nN = len(N)
+    for i in range(nN):
+        for n in N[i]:
+            ofile.write('{}\t'.format(n))
+        ofile.write('\n')
+        for r in rates[i]:
+            ofile.write('{}\t'.format(r))
+        ofile.write('\n')
+    ofile.close()
+
+if __name__=="__main__":
+    L = 25
+    W = 11
+    Nnom = 7
+
+    NAll = gen_N(L,W,Nnom,verbose=True)
+    nN = len(NAll)
+    print(nN)
+
+    Nout = [[] for i in range(nN)]
+    rates = [[] for i in range(nN)]
+
+    def mapFun(N):
+        rates = calc_rates(L,W,N,verbose=True)
+        return rates
+
+    nproc = 30
+    #pool = Pool(processes=nproc)
+    pool = Pool()
+    for ind, res in enumerate(pool.imap(mapFun, NAll, chunksize=3)):
+        Nout[ind] = res[0]
+        rates[ind] = res[1]
+
+    fname = 'L{}W{}.rates'.format(L,W)
+    saveRates(NAll,rates,fname)

setup.py

@@ -36,7 +36,7 @@ def get_version():
     keywords=['fracture mechanics',
               'statistical mechanics',
               'thermodynamics'],
-    install_requires=['matplotlib', 'numpy', 'scipy'],
+    install_requires=['matplotlib', 'numpy', 'scipy', 'sympy'],
     extras_require={
       'docs': ['docutils>=0.14,<0.21', 'ipykernel', 'ipython',
                'jinja2>=3.0', 'nbsphinx', 'pycodestyle',

statmechcrack/mechanical.py

@@ -112,8 +112,8 @@ def beta_U_0(self, v, s):
 
         .. math::
             \beta U_0(v, \mathbf{s}) =
-            \sum_{i=1}^{L-1} \frac{\kappa}{2} \left(
-                s_{i-1} - 2s_i + s_{i+1}
+            \sum_{i=2}^{L} \frac{\kappa}{2} \left(
+                s_{i-2} - 2s_{i-1} + s_i
             \right)^2,
 
         where :math:`s_0\equiv v`.
@@ -746,7 +746,7 @@ def minimize_beta_U_00(self, v, lambda_):
                 - (*numpy.ndarray*) -
                   The minimized nondimensional potential energy.
                 - (*numpy.ndarray*) -
-                  The corresponding nondimensional end separation.
+                  The corresponding nondimensional end force.
                 - (*numpy.ndarray*) -
                   The corresponding nondimensional positions.
 
@@ -807,6 +807,8 @@ def minimize_beta_U(self, v, transition_state=False):
 
                 - (*numpy.ndarray*) -
                   The minimized nondimensional potential energy.
+                - (*numpy.ndarray*) -
+                  The corresponding nondimensional force.
                 - (*numpy.ndarray*) -
                   The corresponding nondimensional positions.
 

statmechcrack/q2d/isometric.py

@@ -3,6 +3,9 @@
 
 """
 
+import numpy as np
+from numpy.linalg import det, inv
+
 from .mechanical import CrackQ2DMechanical
 
 
@@ -19,3 +22,35 @@ def __init__(self, **kwargs):
 
         """
         CrackQ2DMechanical.__init__(self, **kwargs)
+
+    def k_isometric(self, v, transition_state):
+        r"""The nondimensional forward reaction rate coefficient
+        as a function of the nondimensional end separations
+        in the isometric ensemble.
+
+        Args:
+            v (array_like): The nondimensional end separations.
+            transition_state (int): The kth transition state.
+
+        Returns:
+            numpy.ndarray: The nondimensional forward reaction rate.
+
+        """
+        v_ref = np.ones(self.W)
+        beta_U, _, _, hess = self.minimize_beta_U(v)
+        beta_U_ref, _, _, hess_ref = self.minimize_beta_U(v_ref)
+        beta_U_TS, _, _, hess_TS = self.minimize_beta_U(
+            v, transition_state=transition_state
+        )
+        beta_U_TS_ref, _, _, hess_TS_ref = self.minimize_beta_U(
+            v_ref, transition_state=transition_state
+        )
+        scale = self.varepsilon*(self.W*self.L)**(1/3)
+        return np.exp(
+            beta_U - beta_U_ref - beta_U_TS + beta_U_TS_ref
+        )*np.sqrt(
+            det(hess/scale) /
+            det(hess_ref/scale) *
+            det(hess_TS_ref/scale) /
+            det(hess_TS/scale)
+        )

statmechcrack/q2d/isotensional.py

@@ -3,6 +3,9 @@
 
 """
 
+import numpy as np
+from numpy.linalg import det, inv
+
 from .isometric import CrackQ2DIsometric
 
 
@@ -19,3 +22,35 @@ def __init__(self, **kwargs):
 
         """
         CrackQ2DIsometric.__init__(self, **kwargs)
+
+    def k_isotensional(self, p, transition_state):
+        r"""The nondimensional forward reaction rate coefficient
+        as a function of the nondimensional end forces
+        in the isotensional ensemble.
+
+        Args:
+            p (array_like): The nondimensional end forces.
+            transition_state (int): The kth transition state.
+
+        Returns:
+            numpy.ndarray: The nondimensional forward reaction rate.
+
+        """
+        p_ref = np.zeros(self.W)
+        beta_Pi, _, _, hess = self.minimize_beta_Pi(p)
+        beta_Pi_ref, _, _, hess_ref = self.minimize_beta_Pi(p_ref)
+        beta_Pi_TS, _, _, hess_TS = self.minimize_beta_Pi(
+            p, transition_state=transition_state
+        )
+        beta_Pi_TS_ref, _, _, hess_TS_ref = self.minimize_beta_Pi(
+            p_ref, transition_state=transition_state
+        )
+        scale = self.varepsilon*(self.W*self.L)**(1/3)
+        return np.exp(
+            beta_Pi - beta_Pi_ref - beta_Pi_TS + beta_Pi_TS_ref
+        )*np.sqrt(
+            det(hess/scale) /
+            det(hess_ref/scale) *
+            det(hess_TS_ref/scale) /
+            det(hess_TS/scale)
+        )