notebook.community

Edit and run

Introduction

Preperation

Import packages





In [ ]:



    


import numpy as np
import matplotlib.pyplot as plt
import mapp4py
from mapp4py import md
from lib.elasticity import rot, cubic, resize, displace, crack

Block the output of all cores except for one





In [ ]:



    


from mapp4py import mpi
if mpi().rank!=0:
    with open(os.devnull, 'w') as f:
        sys.stdout = f;

Define an md.export_cfg object

md.export_cfg has a call method that we can use to create quick snapshots of our simulation box





In [ ]:



    


xprt = md.export_cfg("");

Asymptotic Displacement Field of Crack from Linear Elasticity





In [ ]:



    


_ = np.array([[-1,1,0],[1,1,1],[1,1,-2]], dtype=np.float);
Q = np.linalg.inv(np.sqrt(_ @ _.T)) @ _;



    











In [ ]:



    


C = rot(cubic(1.3967587463636366,0.787341583191591,0.609615090769241),Q)



    











In [ ]:



    


B = np.linalg.inv(
    np.array([
    [C[0, 0, 0, 0], C[0, 0, 1, 1], C[0, 0, 0, 1]],
    [C[0, 0, 1, 1], C[1, 1, 1, 1], C[1, 1, 0, 1]],
    [C[0, 0, 0, 1], C[1, 1, 0, 1], C[0, 1, 0, 1]]
    ]
))



    











In [ ]:



    


_ = np.roots([B[0, 0], -2.0*B[0, 2],2.0*B[0, 1]+B[2, 2], -2.0*B[1, 2], B[1, 1]])

mu = np.array([_[0],0.0]);

if np.absolute(np.conjugate(mu[0]) - _[1]) > 1.0e-12:
    mu[1] = _[1];
else:
    mu[1] = _[2]

alpha = np.real(mu);
beta = np.imag(mu);

p = B[0,0] * mu**2 - B[0,2] * mu + B[0, 1]
q = B[0,1] * mu - B[0, 2] + B[1, 1]/ mu

K = np.stack([p, q]) * np.array(mu[1], mu[0]) /(mu[1] - mu[0])

K_r = np.real(K)
K_i = np.imag(K)



    











In [ ]:



    


Tr = np.stack([
    np.array(np.array([[1.0, alpha[0]], [0.0, beta[0]]])), 
    np.array([[1.0, alpha[1]], [0.0, beta[1]]])
], axis=1)


def u_f0(x): return np.sqrt(np.sqrt(x[0] * x[0] + x[1] * x[1]) + x[0])
def u_f1(x): return np.sqrt(np.sqrt(x[0] * x[0] + x[1] * x[1]) - x[0]) * np.sign(x[1]) 


def disp(x): 
    _ = Tr @ x
    return K_r @ u_f0(_) + K_i @ u_f1(_)



    











In [ ]:



    


n = 300;
r = 10;
disp_scale = 0.3;

n0 = int(np.round(n/ (1 +np.pi), ))
n1 = n - n0

xs = np.concatenate((
        np.stack([np.linspace(0, -r , n0), np.full((n0,), -1.e-8)]),
        r * np.stack([np.cos(np.linspace(-np.pi, np.pi , n1)),np.sin(np.linspace(-np.pi, np.pi , n1))]), 
        np.stack([np.linspace(-r, 0 , n0), np.full((n0,), 1.e-8)]),
    ), axis =1)

xs_def = xs + disp_scale * disp(xs)

fig, ax = plt.subplots(figsize=(10.5,5), ncols = 2)
ax[0].plot(xs[0], xs[1], "b-", label="non-deformed");
ax[1].plot(xs_def[0], xs_def[1], "r-.", label="deformed");

Configuration

Create a $[\bar{1}10]\times\frac{1}{2}[111]\times[11\bar{2}]$ cell

start with a $[100]\times[010]\times[001]$ cell





In [ ]:



    


sim = md.atoms.import_cfg("configs/Fe_300K.cfg");
a = sim.H[0][0]

Create a $[\bar{1}10]\times[111]\times[11\bar{2}]$ cell





In [ ]:



    


sim.cell_change([[-1,1,0],[1,1,1],[1,1,-2]])

Remove half of the atoms and readjust the position of remaining

Now one needs to cut the cell in half in $[111]$ direction. We can achive this in three steps:

  1. Remove the atoms that are above located above $\frac{1}{2}[111]$
  2. Double the position of the remiaing atoms in the said direction
  3. Shrink the box affinly to half on that direction




In [ ]:



    


H = np.array(sim.H);
def _(x):
    if x[1] > 0.5*H[1, 1] - 1.0e-8:
        return False;
    else:
        x[1] *= 2.0;
sim.do(_);

_ = np.full((3,3), 0.0)
_[1, 1] = -0.5
sim.strain(_)

Readjust the postions





In [ ]:



    


H = np.array(sim.H);
displace(sim,np.array([sim.H[0][0]/6.0, sim.H[1][1]/6.0, sim.H[2][2]/6.0]))

Replicating the unit cell





In [ ]:



    


max_natms=100000
H=np.array(sim.H);
n_per_area=sim.natms/(H[0,0] * H[1,1]);
_ =np.sqrt(max_natms/n_per_area);
N0 = np.array([
    np.around(_ / sim.H[0][0]),
    np.around(_ / sim.H[1][1]), 
    1], dtype=np.int32)

# make sure in 1 direction it is an even number
if N0[1] % 2 == 1:
    N0[1] += 1

sim *= N0;

Add vacuum





In [ ]:



    


vaccum = 100.0



    











In [ ]:



    


H = np.array(sim.H);
H_new = np.array(sim.H);
H_new[0][0] += vaccum
H_new[1][1] += vaccum
resize(sim, H_new, H.sum(axis=0) * 0.5)

Get the displacement field for this configuration





In [ ]:



    


_ = np.array([[-1,1,0],[1,1,1],[1,1,-2]], dtype=np.float);
Q = np.linalg.inv(np.sqrt(_ @ _.T)) @ _;
C = rot(cubic(1.3967587463636366,0.787341583191591,0.609615090769241),Q)
disp = crack(C)

Impose the diplacement field and other boundary conditions





In [ ]:



    


fixed_layer_thickness = 20.0
intensity = 0.5
rate = 0.001



    











In [ ]:



    


H = np.array(sim.H);
ctr = H.sum(axis=0) * 0.5

lim = np.array([H[0, 0], H[1, 1]]) 
lim -= vaccum;
lim *= 0.5
lim -= fixed_layer_thickness

def _(x, x_d, x_dof):
    
    x_rel = x[:2] - ctr[:2]
    u = disp(x_rel)
    x[:2] += intensity * u
    
    if (np.abs(x_rel) < lim).sum() != 2:
        x_d[:2] = rate * u
        x_dof[0] = False;
        x_dof[1] = False;
    
sim.do(_)



    











In [ ]:



    


md.export_cfg("", extra_vecs=["x_dof"] )(sim, "dumps/crack.cfg")

assign intial velocities





In [ ]:



    


sim.kB = 8.617330350e-5
sim.hP = 4.13566766225 * 0.1 * np.sqrt(1.60217656535/1.66053904020)
sim.create_temp(300.0, 846244)

add hydrogen to the system





In [ ]:



    


sim.add_elem('H',1.007940)

define ensemble

muvt





In [ ]:



    


# GPa and Kelvin
def mu(p,T):
    return -2.37+0.0011237850013293155*T+0.00004308665175*T*np.log(p)-0.000193889932875*T*np.log(T);



    











In [ ]:



    


muvt = md.muvt(mu(1.0e-3,300.0), 300.0, 0.1, 'H', 73108204);
muvt.nevery = 100;
muvt.nattempts=40000;
muvt.ntally=1000;
muvt.export=md.export_cfg('dumps/dump',10000)

run gcmc





In [ ]:



    


muvt.run(sim,100000);

Content source: sinamoeini/mapp4py

Similar notebooks:

notebook.community | gallery | about