Custom Symmetry Group Example
import symd
import matplotlib.pyplot as plt
import numpy as np
Make custom group
Rotational group with reflection boundary conditions. Asymmetric unit is just a guess - not actually correct, so we may have overlapping particles and need to re-run the simulation.
rot = 7
genpos = []
for ri in range(rot):
c = np.round(np.cos(ri * np.pi * 2 / rot), 4)
s = np.round(np.sin(ri * np.pi * 2 / rot), 4)
# convert numbers to affine notation
genpos.append(
f"{c}x {-s:+}y {-0.5 * c + 0.5 * s + 0.5:+},{s}x {c:+}y {-0.5 * s - 0.5 * c + 0.5:+}"
)
specpos = [symd.Group("Oblique", ["1/2, 1/2"], specpos=[])]
asymm = f"1/2≤x≤1.0;1/2≤y≤1.0"
my_group = symd.Group("Oblique", genpos, specpos, "quasi", asymm)
print(my_group)
Group(lattice='Oblique', genpos=['1.0x -0.0y +0.0,0.0x +1.0y +0.0', '0.6235x -0.7818y +0.57915,0.7818x +0.6235y -0.20265', '-0.2225x -0.9749y +1.0987,0.9749x -0.2225y +0.12380000000000002', '-0.901x -0.4339y +1.16745,0.4339x -0.901y +0.73355', '-0.901x +0.4339y +0.73355,-0.4339x -0.901y +1.16745', '-0.2225x +0.9749y +0.12380000000000002,-0.9749x -0.2225y +1.0987', '0.6235x +0.7818y -0.20265,-0.7818x +0.6235y +0.57915'], specpos=[Group(lattice='Oblique', genpos=['1/2, 1/2'], specpos=[], name=None, asymm_unit=None)], name='quasi', asymm_unit='1/2≤x≤1.0;1/2≤y≤1.0')
Run
Run the simulation
def run_sim(
n, number_density, group, images, w=None, retries=50, pos_frames=0, steps=10**6
):
for i in range(retries):
print("Trying on ", i)
try:
np.random.seed(i)
cell = symd.groups.get_cell(number_density, group, 2, n, w)
md = symd.Symd(
nparticles=n,
cell=cell,
ndims=2,
images=images,
force="lj",
wyckoffs=w,
group=group,
steps=steps,
exeDir="quasi",
start_temperature=0.5,
temperature=0.1,
pressure=0.25,
)
md.remove_overlap()
md.runParams["box_update_period"] = 10
md.runParams["langevin_gamma"] = 0.5
md.log_positions(frames=pos_frames)
try:
md.run()
except RuntimeError as e:
d = md.number_density()
if d < 0.6:
print("Not dense enough, retrying", d)
continue
# Basically E-min
md.runParams["start_temperature"] = 0.1
md.runParams["temperature"] = 1e-2
md.runParams["langevin_gamma"] = 0.5
md.runParams["Pressure"] = None
md.runParams["box_update_period"] = 0
md.runParams["steps"] = steps // 10
if pos_frames > 0:
md.log_positions(filename="equil.xyz", frames=pos_frames // 10)
try:
md.run()
except RuntimeError as e:
continue
config = md.positions[-1]
break
except RuntimeError as e:
print(e)
md = None
return md
md = run_sim(256, 0.1, my_group, [0, 0], pos_frames=1000, w=[1])
Trying on 0
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[4], line 1
----> 1 md = run_sim(256, 0.1, my_group, [0, 0], pos_frames=1000, w=[1])
Cell In[3], line 8, in run_sim(n, number_density, group, images, w, retries, pos_frames, steps)
6 try:
7 np.random.seed(i)
----> 8 cell = symd.groups.get_cell(number_density, group, 2, n, w)
9 md = symd.Symd(
10 nparticles=n,
11 cell=cell,
(...)
21 pressure=0.25,
22 )
23 md.remove_overlap()
File /opt/hostedtoolcache/Python/3.12.4/x64/lib/python3.12/site-packages/symd/groups.py:653, in get_cell(number_density, group, dim, n, w)
636 def get_cell(
637 number_density: float,
638 group: Union[int, Group],
(...)
641 w: Optional[List[int]] = None,
642 ) -> List[float]:
643 """
644 Compute unit cell given number density, group, and number of particles in asymmetric unit.
645
(...)
651 :return: flattened unit cell (for use in symd MD engine)
652 """
--> 653 import scipy.optimize as opt
655 if w is None:
656 w = []
ModuleNotFoundError: No module named 'scipy'
plt.plot(md.positions[-1, :, 0], md.positions[-1, :, 1], ".")