diff --git a/pyxtal/__init__.py b/pyxtal/__init__.py index bd51dbb9..cf09a6a9 100644 --- a/pyxtal/__init__.py +++ b/pyxtal/__init__.py @@ -1638,7 +1638,6 @@ def build(self, group, species, numIons, lattice, sites, tol=1e-2, dim=3, use_ha # msg = 'The input lattice needs to be a pyxtal.lattice.Lattice class' # raise ValueError(msg, lattice) self.lattice = lattice - self.dim = dim self.factor = 1.0 self.PBC = self.group.PBC @@ -1667,12 +1666,15 @@ def build(self, group, species, numIons, lattice, sites, tol=1e-2, dim=3, use_ha elif len(wp) == 4: # tuple: (key, x, y, z) = wp _wp = choose_wyckoff(self.group, site=key, dim=dim) + #print('debug build', key, x, y, z, _wp.get_label()) if _wp is not False: if _wp.get_dof() == 0: # fixed pos pt = [0.0, 0.0, 0.0] else: ans = _wp.get_all_positions([x, y, z]) pt = ans[0] if ans is not None else None + #print('debug build', ans, x, y, z) + # print('debug', ans) if pt is not None: _sites.append(atom_site(_wp, pt, sp)) else: @@ -3843,22 +3845,27 @@ def from_tabular_representation( # Conversion from discrete to continuous if discrete: + #print('discrete', x, y, z) [x, y, z] = wp.from_discrete_grid([x, y, z], N_grids) + #print('conversion', x, y, z) # ; print(wp.get_label(), xyz) - xyz = wp.search_generator([x, y, z], tol=tol) + xyz = wp.search_generator([x, y, z], tol=tol, symmetrize=True) if xyz is not None: xyz, wp, _ = wp.merge(xyz, np.eye(3), 1e-3) label = wp.get_label() # ; print(x, y, z, label, xyz[0], xyz[1], xyz[2]) sites.append((label, xyz[0], xyz[1], xyz[2])) numIons += wp.multiplicity + if verbose: + print('add wp', label, xyz) else: if verbose: - print("Cannot find generator from", x, y, z) - print(wp) + print("Cannot get wp in", x, y, z, wp.get_label()) + if len(sites) > 0: try: + # print(sites) self.build(group, ["C"], [numIons], lattice, [sites]) except: print("Invalid Build", number, lattice, numIons, sites) diff --git a/pyxtal/database/HM_Full.csv b/pyxtal/database/HM_Full.csv index a356d357..08032221 100644 --- a/pyxtal/database/HM_Full.csv +++ b/pyxtal/database/HM_Full.csv @@ -50,7 +50,7 @@ Hall,Spg_num,Spg_full,Symbol,P,P^-1,Permutation 49,9,9:-c2,A 1 1 n,"-a-c,a,-b","b,-c,-a-b",1 50,9,9:-c3,I 1 1 a,"c,-a-c,-b","-a-b,-c,a",1 51,9,9:a1,B b 1 1,"b,c,a","c,a,b",1 -52,9,9:a2,C n 1 1,"-b,a,c","b,-a,c",1 +52,9,9:a2,C n 1 1,"b,a,-a-c","b,a,-b-c",1 53,9,9:a3,I c 1 1,"b,-a-c,c","-b-c,a,c",1 54,9,9:-a1,C c 1 1,"-b,a,c","b,-a,c",1 55,9,9:-a2,B n 1 1,"-b,-a-c,a","c,-a,-b-c",1 diff --git a/pyxtal/lego/SO3.py b/pyxtal/lego/SO3.py index 9f29b64a..f8f1ff03 100644 --- a/pyxtal/lego/SO3.py +++ b/pyxtal/lego/SO3.py @@ -35,7 +35,8 @@ def __init__(self, nmax=3, lmax=3, rcut=3.5, alpha=2.0, self.norm = np.sqrt(2*np.sqrt(2)*np.pi/np.sqrt(2*self.ls+1)) self.keys = ['keys', '_nmax', '_lmax', '_rcut', '_alpha', '_cutoff_function', 'weight_on', 'neighborcalc', - 'ncoefs', 'ls', 'norm', 'tril_indices'] + 'ncoefs', 'ls', 'norm', 'tril_indices', 'args'] + self.args = (self.nmax, self.lmax, self.rcut, self.alpha, self._cutoff_function) def __str__(self): s = "SO3 descriptor with Cutoff: {:6.3f}".format(self.rcut) @@ -160,20 +161,21 @@ def compute_p(self, atoms, atom_ids=None): p array (N, M) """ + if atom_ids is None: atom_ids = range(len(atoms)) self.init_atoms(atoms, atom_ids) - plist = np.zeros((len(atoms), self.ncoefs), dtype=np.float64) + plist = np.zeros((len(atom_ids), self.ncoefs), dtype=np.float64) if len(self.neighborlist) > 0: - cs = compute_cs(self.neighborlist, self.nmax, self.lmax, self.rcut, self.alpha, self._cutoff_function) + cs = compute_cs(self.neighborlist, *self.args) cs *= self.atomic_weights[:, np.newaxis, np.newaxis, np.newaxis] cs = np.einsum('inlm,l->inlm', cs, self.norm) # Get r_ij and compute C*np.conj(C) - for i in range(len(atoms)): - centers = self.neighbor_indices[:,0] == i + for _i, i in enumerate(atom_ids): + centers = self.neighbor_indices[:, 0] == i if len(centers) > 0: ctot = cs[centers].sum(axis=0) P = np.einsum('ijk,ljk->ilj', ctot, np.conj(ctot)).real - plist[i] = P[self.tril_indices].flatten() + plist[_i] = P[self.tril_indices].flatten() return plist def compute_dpdr(self, atoms, atom_ids=None): @@ -188,13 +190,14 @@ def compute_dpdr(self, atoms, atom_ids=None): dpdr array (N, N, M, 3) and p array (N, M) """ + if atom_ids is None: atom_ids = range(len(atoms)) self.init_atoms(atoms, atom_ids) - p_list = np.zeros((self.natoms, self.ncoefs), dtype=np.float64) - dp_list = np.zeros((self.natoms, self.natoms, self.ncoefs, 3), dtype=np.float64) + p_list = np.zeros((len(atom_ids), self.ncoefs), dtype=np.float64) + dp_list = np.zeros((len(atom_ids), self.natoms, self.ncoefs, 3), dtype=np.float64) if len(self.neighborlist) > 0: # get expansion coefficients and derivatives - cs, dcs = compute_dcs(self.neighborlist, self.nmax, self.lmax, self.rcut, self.alpha, self._cutoff_function) + cs, dcs = compute_dcs(self.neighborlist, *self.args) # weight cs and dcs cs *= self.atomic_weights[:, np.newaxis, np.newaxis, np.newaxis] @@ -203,12 +206,11 @@ def compute_dpdr(self, atoms, atom_ids=None): dcs = np.einsum('inlmj,l->inlmj', dcs, self.norm) #print('cs, dcs', self.neighbor_indices, cs.shape, dcs.shape) - # Assign cs and dcs to P and dP # cs: (N_ij, n, l, m) => P (N_i, N_des) # dcs: (N_ij, n, l, m, 3) => dP (N_i, N_j, N_des, 3) # (n, l, m) needs to be merged to 1 dimension - for i in range(len(atoms)): + for _i, i in enumerate(atom_ids): # find atoms for which i is the center centers = self.neighbor_indices[:, 0] == i @@ -219,7 +221,7 @@ def compute_dpdr(self, atoms, atom_ids=None): # power spectrum P = c*c_conj # eq_3 (n, n', l) eliminate m P = np.einsum('ijk, ljk->ilj', ctot, np.conj(ctot)).real - p_list[i] += P[self.tril_indices].flatten() + p_list[_i] += P[self.tril_indices].flatten() # gradient of P for each neighbor, eq_26 # (N_ijs, n, n', l, 3) @@ -233,12 +235,13 @@ def compute_dpdr(self, atoms, atom_ids=None): # QZ: to check ijs = self.neighbor_indices[centers] for _id, j in enumerate(ijs[:, 1]): - dp_list[i, j, :, :] += dP[_id][self.tril_indices].flatten().reshape(self.ncoefs, 3) - dp_list[i, i, :, :] -= dP[_id][self.tril_indices].flatten().reshape(self.ncoefs, 3) + tmp = dP[_id][self.tril_indices].flatten().reshape(self.ncoefs, 3) + dp_list[_i, j, :, :] += tmp + dp_list[_i, i, :, :] -= tmp return dp_list, p_list - def compute_dpdr_5d(self, atoms): + def compute_dpdr_5d(self, atoms, atom_ids=None): """ Compute the powerspectrum function with respect to supercell @@ -248,57 +251,52 @@ def compute_dpdr_5d(self, atoms): Returns: dpdr array (N, N, M, 3, 27) and p array (N, M) """ + if atom_ids is None: atom_ids = range(len(atoms)) + self.init_atoms(atoms, atom_ids) + p_list = np.zeros((len(atom_ids), self.ncoefs), dtype=np.float64) + dp_list = np.zeros((len(atom_ids), self.natoms, self.ncoefs, 3, 27), dtype=np.float64) - self.init_atoms(atoms) - p_list = np.zeros((self.natoms, self.ncoefs), dtype=np.float64) - dp_list = np.zeros((self.natoms, self.natoms, self.ncoefs, 3, 27), dtype=np.float64) - - # get expansion coefficients and derivatives - cs, dcs = compute_dcs(self.neighborlist, self.nmax, self.lmax, self.rcut, self.alpha, self._cutoff_function) - - # weight cs and dcs - cs *= self.atomic_weights[:, np.newaxis, np.newaxis, np.newaxis] - dcs *= self.atomic_weights[:, np.newaxis, np.newaxis, np.newaxis, np.newaxis] - cs = np.einsum('inlm,l->inlm', cs, self.norm) - dcs = np.einsum('inlmj,l->inlmj', dcs, self.norm) - #print('cs, dcs', self.neighbor_indices, cs.shape, dcs.shape) - - # Assign cs and dcs to P and dP - # cs: (N_ij, n, l, m) => P (N_i, N_des) - # dcs: (N_ij, n, l, m, 3) => dP (N_i, N_j, N_des, 3) - # (n, l, m) needs to be merged to 1 dimension - neigh_ids = np.arange(len(self.neighbor_indices)) - for i in range(len(atoms)): - # find atoms for which i is the center - pair_ids = neigh_ids[self.neighbor_indices[:, 0] == i] - if len(pair_ids) > 0: - ctot = cs[pair_ids].sum(axis=0) #(n, l, m) - dctot = dcs[pair_ids].sum(axis=0) - # power spectrum P = c*c_conj - # eq_3 (n, n', l) eliminate m - P = np.einsum('ijk, ljk->ilj', ctot, np.conj(ctot)).real - p_list[i] += P[self.tril_indices].flatten() - - # loop over each pair - for pair_id in pair_ids: - (_, j, x, y, z) = self.neighbor_indices[pair_id] - # map from (x, y, z) to (0, 27) - cell_id = (x+1) * 9 + (y+1) * 3 + z + 1 + if len(self.neighborlist) > 0: + # get expansion coefficients and derivatives + cs, dcs = compute_dcs(self.neighborlist, *self.args) - # gradient of P for each neighbor, eq_26 - # (N_ijs, n, n', l, 3) - # dc * c_conj + c * dc_conj - dP = np.einsum('ijkn, ljk->iljn', dcs[pair_id], np.conj(ctot)) - dP += np.einsum('ijkn, ljk->iljn', np.conj(dcs[pair_id]), ctot) - #dP += np.conj(np.transpose(dP, axes=[1, 0, 2, 3])) - #dP += np.einsum('ijkn, ljk->iljn', np.conj(dctot), cs[pair_id]) - #dP += np.einsum('ijkn, ljk->iljn', dctot, np.conj(cs[pair_id])) + # weight cs and dcs + cs *= self.atomic_weights[:, np.newaxis, np.newaxis, np.newaxis] + dcs *= self.atomic_weights[:, np.newaxis, np.newaxis, np.newaxis, np.newaxis] + cs = np.einsum('inlm,l->inlm', cs, self.norm) + dcs = np.einsum('inlmj,l->inlmj', dcs, self.norm) + #print('cs, dcs', self.neighbor_indices, cs.shape, dcs.shape) + + # cs: (N_ij, n, l, m) => P (N_i, N_des) + # dcs: (N_ij, n, l, m, 3) => dP (N_i, N_j, N_des, 3) + # (n, l, m) needs to be merged to 1 dimension + neigh_ids = np.arange(len(self.neighbor_indices)) + for _i, i in enumerate(atom_ids): + # find atoms for which i is the center + pair_ids = neigh_ids[self.neighbor_indices[:, 0] == i] + if len(pair_ids) > 0: + ctot = cs[pair_ids].sum(axis=0) #(n, l, m) + dctot = dcs[pair_ids].sum(axis=0) + # power spectrum P = c * c_conj + P = np.einsum('ijk, ljk->ilj', ctot, np.conj(ctot)).real + p_list[_i] += P[self.tril_indices].flatten() + + # loop over each pair + for pair_id in pair_ids: + (_, j, x, y, z) = self.neighbor_indices[pair_id] + # map from (x, y, z) to (0, 27) + cell_id = (x+1) * 9 + (y+1) * 3 + z + 1 + + # (N_ijs, n, n', l, 3) + # dp = dc * c_conj + c * dc_conj + dP = np.einsum('ijkn, ljk->iljn', dcs[pair_id], np.conj(ctot)) + dP += np.einsum('ijkn, ljk->iljn', np.conj(dcs[pair_id]), ctot) - dP = dP.real[self.tril_indices].flatten().reshape(self.ncoefs, 3) - #print(cs[pair_id].shape, dcs[pair_id].shape, dP.shape) + dP = dP.real[self.tril_indices].flatten().reshape(self.ncoefs, 3) + #print(cs[pair_id].shape, dcs[pair_id].shape, dP.shape) - dp_list[i, j, :, :, cell_id] += dP - dp_list[i, i, :, :, 13] -= dP + dp_list[_i, j, :, :, cell_id] += dP + dp_list[_i, i, :, :, 13] -= dP return dp_list, p_list diff --git a/pyxtal/lego/builder.py b/pyxtal/lego/builder.py index c28af612..ca0d885f 100644 --- a/pyxtal/lego/builder.py +++ b/pyxtal/lego/builder.py @@ -39,6 +39,7 @@ from time import time # np.set_printoptions(precision=3, suppress=True) +VECTORS = np.array([[x1, y1, z1] for x1 in range(-1, 2) for y1 in range(-1, 2) for z1 in range(-1, 2)]) def generate_wp_lib_par(spgs, composition, num_wp, num_fu, num_dof): """ @@ -53,7 +54,6 @@ def generate_wp_lib_par(spgs, composition, num_wp, num_fu, num_dof): my_spgs.append(spg) return (my_spgs, wp_libs) - def generate_xtal_par(wp_libs, niter, dim, elements, calculator, ref_environments, criteria, T, N_max, early_quit): """ @@ -70,7 +70,6 @@ def generate_xtal_par(wp_libs, niter, dim, elements, calculator, ref_environment return (xtals, sims) - def minimize_from_x_par(*args): """ A wrapper to call minimize_from_x function in parallel @@ -141,7 +140,9 @@ def generate_xtal(dim, spg, wps, niter, elements, calculator, def minimize_from_x(x, dim, spg, wps, elements, calculator, ref_environments, T=0.2, niter=20, early_quit=0.02, opt_type='local', minimizers=[('Nelder-Mead', 100), ('L-BFGS-B', 100)], - filename='local_opt_data.txt', random_state=None): + filename='local_opt_data.txt', random_state=None, + #derivative=True): + derivative=False): """ Generate xtal from the 1d representation @@ -151,6 +152,11 @@ def minimize_from_x(x, dim, spg, wps, elements, calculator, ref_environments, wps (string): e.g. [['4a', '8b']] elements (string): e.g., ['Si', 'O'] """ + if derivative: + jac = calculate_dSdx + else: + jac = None + g, wps, dof = get_input_from_letters(spg, wps, dim) l_type = g.lattice_type sites_wp = [] @@ -208,7 +214,7 @@ def minimize_from_x(x, dim, spg, wps, elements, calculator, ref_environments, except: return None - x0 = x.copy() + x0 = np.array(x.copy()) # Extract variables, call from Pyxtal [N_abc, N_ang] = Lattice.get_dofs(xtal.lattice.ltype) rep = xtal.get_1D_representation() @@ -278,6 +284,7 @@ def print_local_fun(x): res = minimize(calculate_S, x, method=method, args=(xtal, ref_envs, calculator), + jac=None if method=='Nelder-Mead' else jac, bounds=bounds, options={'maxiter': step}, # 'disp': True}, callback=callback) @@ -308,11 +315,10 @@ def print_fun_local(x): 'fatol': 1e-6, 'ftol': 1e-6}} - bounded_step = RandomDisplacementBounds(np.array([b[0] for b in bounds]), - np.array([b[1] - for b in bounds]), - id1=N_abc + N_ang, - id2=N_abc) + bounded_step = RandomDispBounds(np.array([b[0] for b in bounds]), + np.array([b[1] for b in bounds]), + id1=N_abc + N_ang, + id2=N_abc) # set call back function for debugging def print_fun(x, f, accepted): @@ -358,38 +364,47 @@ def calculate_dSdx(x, xtal, des_ref, f, eps=1e-4, symmetry=True, verbose=False): verbose (bool): output more information """ xtal.update_from_1d_rep(x) - atoms = xtal.to_ase(resort=False, add_vaccum=False) # * 2 - atoms.set_positions(atoms.positions + 1e-3 * - np.random.random([len(atoms), 3])) - ref_pos = atoms.positions - - # results = f.calculate(atoms, ids, derivative=True) - results = f.calculate(atoms, derivative=True) - dPdr = results['dxdr'] # ; print(dpdr>0) - P = results['x'] - - # Compute dSdr - dSdr = np.einsum("ik, ijkl -> jl", 2*(P - des_ref), dPdr) + ids = [0] + weights = [] + for site in xtal.atom_sites: + ids.append(site.wp.multiplicity + ids[-1]) + weights.append(site.wp.multiplicity) + ids = ids[:-1] + weights = np.array(weights, dtype=float) + + atoms = xtal.to_ase() + dPdr, P = f.compute_dpdr_5d(atoms, ids) + dPdr = np.einsum('i, ijklm -> ijklm', weights, dPdr) + + # Compute dSdr [N, M] [N, N, M, 3, 27] => [N, 3, 27] + # print(P.shape, des_ref.shape, dPdr.shape) + dSdr = np.einsum("ik, ijklm -> jlm", 2*(P - des_ref), dPdr) + + # Get supercell positions + ref_pos = np.repeat(atoms.positions[:, :, np.newaxis], 27, axis=2) + cell_shifts = np.dot(VECTORS, atoms.cell) + ref_pos += cell_shifts.T[np.newaxis, :, :] # Compute drdx via numerical func - drdx = np.zeros([len(atoms), 3, len(x)]) + drdx = np.zeros([len(atoms), 3, 27, len(x)]) + xtal0 = xtal.copy() for i in range(len(x)): x0 = x.copy() x0[i] += eps + # Maybe expensive Reduce calls, just need positions xtal0.update_from_1d_rep(x0) - # print(xtal0) - pos = xtal0.to_ase(resort=False, add_vaccum=False).positions - drdx[:, :, i] = (pos - ref_pos)/eps - # print("drdx", i, x0, '\n', drdx[:, :, i]) - - # dPdx = np.einsum('ijkl, klm->ijm', dPdr, drdx) - # dSdx = np.einsum('ij, ijk->k', 2*(P-des_ref), dPdx) - # dSdx = dSdr * drdx - dSdx = np.einsum("ij, ijk -> k", dSdr, drdx) + atoms = xtal0.to_ase() - return dSdx, dSdr, dPdr + # Get supercell positions + pos = np.repeat(atoms.positions[:, :, np.newaxis], 27, axis=2) + cell_shifts = np.dot(VECTORS, atoms.cell) + pos += cell_shifts.T[np.newaxis, :, :] + drdx[:, :, :, i] += (pos - ref_pos)/eps + # [N, 3, 27] [N, 3, 27, H] => H + dSdx = np.einsum("ijk, ijkl -> l", dSdr, drdx) + return dSdx def calculate_S(x, xtal, des_ref, f, verbose=False): """ @@ -415,14 +430,9 @@ def calculate_S(x, xtal, des_ref, f, verbose=False): weights.append(site.wp.multiplicity) ids = ids[:-1] weights = np.array(weights, dtype=float) - # try: - if True: - atoms = xtal.to_ase(resort=False, add_vaccum=False) - des = f.calculate(atoms, ids)['x'][ids] - # except: - # #print("bug in xtal2ase or bad structure", xtal.lattice) - # des = np.zeros(des_ref.shape) - # #raise ValueError('Skip the random error') + atoms = xtal.to_ase(resort=False) + #des = f.calculate(atoms, ids)['x'][ids] + des = f.compute_p(atoms, ids) else: des = np.zeros(des_ref.shape) weights = 1 @@ -618,12 +628,10 @@ def set_descriptor_calculator(self, dtype='SO3', mykwargs={}): 'rcut': 2.2, 'alpha': 1.5, 'weight_on': True, - #'neighborlist': 'ase', } kwargs.update(mykwargs) self.calculator = SO3(**kwargs) - # print(self.calculator)#; import sys; sys.exit() def set_reference_enviroments(self, cif_file, substitute=None): """ @@ -655,12 +663,12 @@ def set_reference_enviroments(self, cif_file, substitute=None): if self.verbose: print("ids from Reference xtal", ids) atoms = xtal.to_ase(resort=False) - self.ref_environments = self.calculator.calculate(atoms, ids)['x'][ids] + self.ref_environments = self.calculator.compute_p(atoms, ids) if self.verbose: print(self.ref_environments) self.ref_xtal = xtal - def set_criteria(self, CN=None, dimension=3, min_density=None, exclude_ii=False): + def set_criteria(self, CN=None, dimension=None, min_density=None, exclude_ii=False): """ define the criteria to check if a structure is good @@ -795,31 +803,12 @@ def optimize_xtals_serial(self, xtals, args): xtals_opt.append(xtal) return xtals_opt - def optimize_reps(self, reps, ncpu=1, opt_type='local', - T=0.2, niter=20, early_quit=0.02, - add_db=True, symmetrize=False, - minimizers=[('Nelder-Mead', 100), ('L-BFGS-B', 100)], - ): - """ - Perform optimization for each structure - - Args: - reps: list of reps - ncpu (int): - """ - args = (opt_type, T, niter, early_quit, add_db, symmetrize, minimizers) - if ncpu > 1: - valid_xtals = self.optimize_reps_mproc(reps, ncpu, args) - return valid_xtals - else: - raise NotImplementedError("optimize_reps works in parallel mode") - - def optimize_reps_mproc(self, reps, ncpu, args): + def optimize_xtals_mproc(self, xtals, ncpu, args): """ Optimization in multiprocess mode. Args: - reps: list of reps + xtals: list of xtals ncpu (int): number of parallel python processes args: (opt_type, T, n_iter, early_quit, add_db, symmetrize, minimizers) """ @@ -830,11 +819,10 @@ def optimize_reps_mproc(self, reps, ncpu, args): # Split the input structures to minibatches N_batches = 10 * ncpu - for _i, i in enumerate(range(0, len(reps), N_batches)): - start, end = i, min([i+N_batches, len(reps)]) + for _i, i in enumerate(range(0, len(xtals), N_batches)): + start, end = i, min([i+N_batches, len(xtals)]) ids = list(range(start, end)) print(f"Rank {self.rank} minibatch {start} {end}") - self.logging.info(f"Rank {self.rank} minibatch {start} {end}") self.print_memory_usage() def generate_args(): @@ -845,12 +833,10 @@ def generate_args(): _ids = ids[j::ncpu] wp_libs = [] for id in _ids: - rep = reps[id] - xtal = pyxtal() - xtal.from_tabular_representation(rep, normalize=False) + xtal = xtals[id] x = xtal.get_1d_rep_x() spg, wps, _ = self.get_input_from_ref_xtal(xtal) - wp_libs.append((x, spg, wps)) + wp_libs.append((x, xtal.group.number, wps)) yield (self.dim, wp_libs, self.elements, self.calculator, self.ref_environments, opt_type, T, niter, early_quit, minimizers) @@ -875,15 +861,59 @@ def generate_args(): gc.collect() # Explicitly call garbage collector to free memory xtals_opt = list(xtals_opt) - print(f"Rank {self.rank} finish optimize_reps_mproc {len(xtals_opt)}") + print(f"Rank {self.rank} finish optimize_xtals_mproc {len(xtals_opt)}") return xtals_opt - def optimize_xtals_mproc(self, xtals, ncpu, args): + def optimize_reps(self, reps, ncpu=1, opt_type='local', + T=0.2, niter=20, early_quit=0.02, + add_db=True, symmetrize=False, + minimizers=[('Nelder-Mead', 100), ('L-BFGS-B', 100)], + discrete=False): + """ + Perform optimization for each structure + + Args: + reps: list of reps + ncpu (int): + """ + args = (opt_type, T, niter, early_quit, add_db, symmetrize, minimizers) + if ncpu == 1: + valid_xtals = self.optimize_reps_serial(reps, args, discrete) + else: + valid_xtals = self.optimize_reps_mproc(reps, ncpu, args, discrete) + return valid_xtals + + def optimize_reps_serial(self, reps, args, discrete): """ Optimization in multiprocess mode. Args: - xtals: list of xtals + reps: list of reps + ncpu (int): number of parallel python processes + args: (opt_type, T, n_iter, early_quit, add_db, symmetrize, minimizers) + """ + xtals_opt = [] + for i, rep in enumerate(reps): + #print('start', i, rep, len(rep)) + xtal = pyxtal() + xtal.from_tabular_representation(rep, + normalize=False, + discrete=discrete) + #print(xtal.get_xtal_string()) + #print(xtal) + xtal, sim, _xs = self.optimize_xtal(xtal, i, *args) + if xtal is not None: + xtals_opt.append(xtal) + else: + print("Debug===="); import sys; sys.exit() + return xtals_opt + + def optimize_reps_mproc(self, reps, ncpu, args, discrete): + """ + Optimization in multiprocess mode. + + Args: + reps: list of reps ncpu (int): number of parallel python processes args: (opt_type, T, n_iter, early_quit, add_db, symmetrize, minimizers) """ @@ -894,10 +924,11 @@ def optimize_xtals_mproc(self, xtals, ncpu, args): # Split the input structures to minibatches N_batches = 10 * ncpu - for _i, i in enumerate(range(0, len(xtals), N_batches)): - start, end = i, min([i+N_batches, len(xtals)]) + for _i, i in enumerate(range(0, len(reps), N_batches)): + start, end = i, min([i+N_batches, len(reps)]) ids = list(range(start, end)) print(f"Rank {self.rank} minibatch {start} {end}") + self.logging.info(f"Rank {self.rank} minibatch {start} {end}") self.print_memory_usage() def generate_args(): @@ -908,10 +939,14 @@ def generate_args(): _ids = ids[j::ncpu] wp_libs = [] for id in _ids: - xtal = xtals[id] + rep = reps[id] + xtal = pyxtal() + xtal.from_tabular_representation(rep, + normalize=False, + discrete=discrete) x = xtal.get_1d_rep_x() spg, wps, _ = self.get_input_from_ref_xtal(xtal) - wp_libs.append((x, xtal.group.number, wps)) + wp_libs.append((x, spg, wps)) yield (self.dim, wp_libs, self.elements, self.calculator, self.ref_environments, opt_type, T, niter, early_quit, minimizers) @@ -936,7 +971,7 @@ def generate_args(): gc.collect() # Explicitly call garbage collector to free memory xtals_opt = list(xtals_opt) - print(f"Rank {self.rank} finish optimize_xtals_mproc {len(xtals_opt)}") + print(f"Rank {self.rank} finish optimize_reps_mproc {len(xtals_opt)}") return xtals_opt def optimize_xtal(self, xtal, count=0, opt_type='local', @@ -1295,13 +1330,8 @@ def process_xtal(self, xtal, sim, count=0, xs=None, energy=None, 'similarity': sim[1], } if xs is not None: - try: - kvp['x_init'] = np.array2string(xs[0]) - kvp['x_opt'] = np.array2string(xs[1]) - except: - print("Error in xs", xs) - kvp['x_init'] = 'N/A' - kvp['x_opt'] = 'N/A' + kvp['x_init'] = np.array2string(xs[0]) + kvp['x_opt'] = np.array2string(xs[1]) if energy is not None: kvp['ff_energy'] = energy if topology is not None: @@ -1315,7 +1345,7 @@ def process_xtal(self, xtal, sim, count=0, xs=None, energy=None, """ -class RandomDisplacementBounds(object): +class RandomDispBounds(object): """ random displacement with bounds: see: https://stackoverflow.com/a/21967888/2320035 diff --git a/pyxtal/symmetry.py b/pyxtal/symmetry.py index 758117ad..cc8c35b4 100644 --- a/pyxtal/symmetry.py +++ b/pyxtal/symmetry.py @@ -2918,7 +2918,7 @@ def search_generator_dist(self, pt, lattice=None, group=None): min_index = np.argmin(distances) return pts[min_index], np.min(distances) - def search_generator(self, pos, ops=None, tol=1e-2): + def search_generator(self, pos, ops=None, tol=1e-2, symmetrize=False): """ search generator for a special Wyckoff position @@ -2926,6 +2926,7 @@ def search_generator(self, pos, ops=None, tol=1e-2): pos: initial xyz position ops: list of symops tol: tolerance + symmetrize (bool): apply symmetrization Return: pos1: the position that matchs the standard setting @@ -2945,6 +2946,8 @@ def search_generator(self, pos, ops=None, tol=1e-2): if diff.sum() < tol: pos1 -= np.floor(pos1) match = True + if symmetrize: + pos1 = pos0 break if match: