add recp codes

Author: qzhu2017

pyxtal/miscellaneous/RECP/dia.py

@@ -0,0 +1,21 @@
+from pyxtal.xtal import RECP
+from pyxtal import pyxtal
+import numpy as np
+
+xtal1 = pyxtal(); xtal1.from_prototype('diamond') #; xtal.to_file('dia.cif'); print(xtal)
+xtal2 = pyxtal(); xtal2.from_prototype('cBN');
+xtal3 = pyxtal(); xtal3.from_prototype('h-diamond') #; xtal.to_file('dia.cif'); print(xtal)
+data = [(xtal1, 'dia'), (xtal2, 'cBN'), (xtal3, 'h-dia')]
+for i in range(3, 5):
+    xtal = pyxtal()
+    xtal.from_seed(f'sub{i}.cif')
+    data.append((xtal, f'sub{i}'))
+
+data0 = []
+recp = RECP(dmax=10.0, nmax=20, lmax=8)
+for i, (xtal, label) in enumerate(data):
+    d, rdf = recp.compute(xtal.to_ase(), norm=True)
+    if i > 0:
+        print(f"{label:10s} : {np.linalg.norm(d - data0[0][0]):12.4f}")
+    data0.append((d, rdf, label))
+recp.plot(data0, filename='reciprocal_dia.png')

pyxtal/miscellaneous/RECP/test0.py

@@ -0,0 +1,54 @@
+from pyxtal.xtal import RECP
+from pyxtal import pyxtal
+import matplotlib.pyplot as plt
+
+# Use a modern style
+plt.style.use('seaborn-v0_8-darkgrid')
+
+# Customize plot appearance
+plt.rcParams.update({
+    #'font.family': 'serif',
+    'font.size': 14,
+    'axes.labelsize': 15,
+    'axes.titlesize': 15,
+    'legend.fontsize': 14,
+    'axes.grid': True,
+    'grid.alpha': 0.3,
+    'lines.linewidth': 2,
+    'figure.dpi': 300
+})
+
+recp = RECP(dmax=10.0, nmax=10, lmax=10, rbasis='Bessel')
+prototypes = ['diamond',
+              'h-diamond',
+              #'graphite',
+              #'$\alpha$-boron',
+              'a-quartz',
+              'b-quartz']
+
+fig, axs = plt.subplots(len(prototypes), 1, figsize=(6, 1.6*len(prototypes)))
+data = []
+for row, prototype in enumerate(prototypes):
+    xtal = pyxtal()
+    xtal.from_prototype(prototype)
+    p, rdf = recp.compute(xtal.to_ase(), norm=True)
+    # Plot p values
+    if prototype == 'a-quartz':
+        label = r'$\alpha$-quartz'
+    elif prototype == 'b-quartz':
+        label = r'$\beta$-quartz'
+    else:
+        label = prototype
+    label += f" ({xtal.group.number}, {xtal.group.symbol})"
+    axs[row].plot(p, label=label)
+    axs[row].legend(loc=1)
+    axs[row].set_ylabel('$P_{nl}$')
+
+    if row == len(prototypes) - 1:
+        axs[row].set_xlabel('Power Spectrum Index')
+    else:
+        axs[row].set_xticklabels([])
+    axs[row].set_ylim(0, 1.0)
+plt.tight_layout()
+plt.savefig('p-demo.png', dpi=300)
+plt.close()

pyxtal/miscellaneous/RECP/test1.py

@@ -0,0 +1,90 @@
+from pyxtal import pyxtal
+from pyxtal.reciprocal import RECP
+import matplotlib.pyplot as plt
+import numpy as np
+np.random.seed(42)
+
+# Use a modern style
+plt.style.use('seaborn-v0_8-darkgrid')
+
+# Customize plot appearance
+plt.rcParams.update({
+    #'font.family': 'serif',
+    'font.size': 12,
+    'axes.labelsize': 14,
+    'axes.titlesize': 14,
+    'legend.fontsize': 12,
+    'axes.grid': True,
+    'grid.alpha': 0.3,
+    'lines.linewidth': 2,
+    'figure.dpi': 300
+})
+
+data = [(227, ['8a'], [3.567]),
+        (216, ['4a', '4d'], [3.567]),
+        (210, ['8b'], [3.567]),
+        (203, ['8a'], [3.567]),
+        (166, ['6c'], [2.522, 6.178, 0.125]),
+        (141, ['4a'], [2.522, 3.567]),
+        (122, ['4b'], [2.522, 3.567]),
+        (119, ['2b', '2d'], [2.522, 3.567]),
+        (109, ['4a'], [2.522, 3.567, 0.875]),
+        (98, ['4b'], [2.522, 3.567]),
+        (88, ['4a'], [2.522, 3.567], ),
+        #(74, ['4e'], [2.522, 2.522, 3.567, 0.875]),
+        (70, ['8a'], [3.567, 3.567, 3.567])]
+
+data0 = []
+recp = RECP(dmax=10.0, nmax=10, lmax=10, rbasis='Bessel')
+
+
+xtal = pyxtal()
+xtal.from_prototype('diamond')
+xtal_sub = xtal.subgroup_once(H=141, eps=0)
+
+fig, axs = plt.subplots(3, 2, figsize=(12, 7.5))
+for row, eps in enumerate([0, 0.02, 0.05]):
+    data1 = []
+    for i, d in enumerate(data):
+        if i == 0:
+            xtal0 = xtal.copy()
+        else:
+            xtal0 = xtal.subgroup_once(H=d[0], eps=eps)
+        if xtal0 is not None:
+            p, rdf = recp.compute(xtal0.to_ase(), norm=True)
+        else:
+            xtal0 = xtal_sub.subgroup_once(H=d[0], eps=eps)
+            if xtal0 is None:
+                print('problem', d[0]); import sys; sys.exit()
+            p, rdf = recp.compute(xtal0.to_ase(), norm=True)
+        data1.append((p, rdf, f'{xtal0.group.number}'))
+
+    # Plot p values
+    axs[row, 0].set_title(f'$P$ (noise={eps})', y=0.80, x=0.05, loc='left')
+    for p, _, label in data1:
+        axs[row, 0].plot(p, label=label, alpha=0.5, lw=1.0)
+
+
+    axs[row, 0].set_ylabel('$P_{nl}$')
+    if row == 2:
+        axs[row, 0].set_xlabel('Power Spectrum Index')
+
+    # Plot rdf values
+    axs[row, 1].set_title(f'$G$ (noise={eps})', y=0.80, x=0.05, loc='left')
+    for _, rdf, label in data1:
+        axs[row, 1].plot(rdf, label=label, alpha=0.5, lw=1.0)
+    if row == 2:
+        axs[row, 1].set_xlabel(r'$d$ ($\mathrm{\AA}^{-1}$)')
+    else:
+        axs[row, 0].set_xticklabels([])
+        axs[row, 1].set_xticklabels([])
+
+    axs[row, 1].set_ylabel('$G(d)$')
+    if row == 0: axs[row, 1].legend(ncol=2, loc='upper right')
+
+plt.tight_layout()
+plt.savefig('test1-diamond.png')
+plt.close()
+
+
+

pyxtal/reciprocal.py

@@ -348,7 +348,7 @@ def apply_bounds(tensor):
                 print(f"Step {step}, {loss_sum:.6f}, LR={scheduler.get_last_lr()[0]:.6f}")
             if step + 1 == num_steps:
                 print(f"stopping at last iteration")
-        xtal =
+        #xtal =
         return rep.detach(), losses.detach()
 
     def loss(self, spg, wps, elements, P_ref, RDF_ref):