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
In [ ]:
from mapp4py import mpi
if mpi().rank!=0:
with open(os.devnull, 'w') as f:
sys.stdout = f;
md.export_cfg
has a call method that we can use to create quick snapshots of our simulation box
In [ ]:
xprt = md.export_cfg("");
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");
In [ ]:
sim = md.atoms.import_cfg("configs/Fe_300K.cfg");
a = sim.H[0][0]
In [ ]:
sim.cell_change([[-1,1,0],[1,1,1],[1,1,-2]])
Now one needs to cut the cell in half in $[111]$ direction. We can achive this in three steps:
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(_)
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]))
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;
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)
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)
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")
In [ ]:
sim.kB = 8.617330350e-5
sim.hP = 4.13566766225 * 0.1 * np.sqrt(1.60217656535/1.66053904020)
sim.create_temp(300.0, 846244)
In [ ]:
sim.add_elem('H',1.007940)
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)
In [ ]:
muvt.run(sim,100000);