diff --git a/pyxtal/__init__.py b/pyxtal/__init__.py index f629c6d9..2a7bb97b 100644 --- a/pyxtal/__init__.py +++ b/pyxtal/__init__.py @@ -1085,12 +1085,12 @@ def _get_formula(self): """ A quick function to get the formula. """ + from pyxtal.database.element import Element formula = "" if self.molecular: numspecies = self.numMols species = [str(mol) for mol in self.molecules] - #print(species, numspecies) else: specie_list = [] for site in self.atom_sites: @@ -1107,6 +1107,7 @@ def _get_formula(self): self.numIons = numIons numspecies = self.numIons for i, s in zip(numspecies, species): + if type(s) == int: s = Element(s).short_name formula += "{:s}{:d}".format(s, int(i)) self.formula = formula @@ -2974,3 +2975,83 @@ def to_atomic_xtal(self): return xtal else: raise RuntimeError('The input must be molecular xtal') + + + def get_forcefield(self, ff_style='openff', code='lammps', + chargemethod='am1bcc'): + """ + An interface to create forcefield for molecular simulation with + - Charmm + - LAMMPS + Note that the current approach generates charge with its own conformer, + need to provide the option to use the conformer from the given molecule + + Args: + - ff_style: 'gaff' or 'openff' + - code: 'lammps' or 'charmm' + - charge_method: 'am1bcc', 'am1-mulliken', 'mmff94', 'gasteiger' + + Returns: + An ase_atoms_objects with force field information + """ + from ost.forcefield import forcefield + from ost.lmp import LAMMPSStructure + from ost.charmm import CHARMMStructure + from ost.interfaces.parmed import ParmEdStructure + + if self.molecular: + smiles = [mol.smile for mol in self.molecules] + assert(smiles[0] is not None) + + # Lammps needs the xyz for all molecules + if code == 'lammps': + # reset lammps cell + from pyxtal.util import reset_lammps_cell + atoms = self.to_ase(resort=False) + atoms = reset_lammps_cell(atoms) + converter = LAMMPSStructure + n_mols = self.numMols + + # Charmm only needs the xyz for the asymmetric unit + else: + converter = CHARMMStructure + coords, species = [], [] + checked = [] + for site in self.mol_sites: + if site.type not in checked: + _coords, _species = site._get_coords_and_species( + absolute=True, first=True) + coords.extend(_coords) + species.extend(_species) + checked.append(site.type) + n_mols = [1] * len(self.molecules) + atoms = Atoms(symbols=species, positions=coords) + + # print(len(atoms), n_mols) + # Initialize the forcefield instance from ost + ff = forcefield(smiles, style=ff_style, chargemethod=chargemethod) + + # Create the Parmed to handle FF parameters + if sum(n_mols) == 1: + pd_struc = ff.molecules[0].copy(cls=ParmEdStructure) + else: + from functools import reduce + from operator import add + mols = [] + for i, m in enumerate(n_mols): + mols += [ff.molecules[i] * m] + pd_struc = reduce(add, mols) + pd_struc.update(atoms) + + ase_with_ff = converter.from_structure(pd_struc) + if code == 'lammps': + print("\n Write lmp.in and lmp.dat ......") + ase_with_ff.write_lammps() + else: + print("\n Write charmm.rtf and charmm.prm ......") + ase_with_ff.write_charmmfiles() + #ase_with_ff.write_charmmfiles_by_parmd() + + return ase_with_ff + else: + raise NotImplementedError("\nNo support for atomic crystals") diff --git a/pyxtal/io.py b/pyxtal/io.py index 5d08c03a..2dc41dc2 100644 --- a/pyxtal/io.py +++ b/pyxtal/io.py @@ -382,7 +382,7 @@ def resort(self, molecules): self.wps.append(wps[i]) self.numMols[j] += len(wps[i]) #print("================================================ADDDDDD", id, len(mol1)) - + # check if some molecules cannot be matched if len(ids_done) < len(ids): #print("==========================================================Nonmatched molecules", ids_done, ids) @@ -431,7 +431,7 @@ def add_Hydrogens(self, smile, xyz): m3 = AllChem.ConstrainedEmbed(m1, m2) except ValueError as e: raise ReadSeedError(str(e)) - + conf = m3.GetConformer(0) #; print(conf.GetPositions()) return conf.GetPositions() diff --git a/pyxtal/molecule.py b/pyxtal/molecule.py index 9ffbfd91..7f07354c 100644 --- a/pyxtal/molecule.py +++ b/pyxtal/molecule.py @@ -1084,6 +1084,12 @@ def get_orientation(self, xyz, rtol=0.15): ang = r.as_euler('zxy', degrees=True) return ang, 0, False + def to_ase(self): + """ + Convert to ase atoms + """ + return self.mol.to_ase_atoms() + def reset_positions(self, coors): """ reset the coordinates diff --git a/pyxtal/util.py b/pyxtal/util.py index 482258eb..69ec5182 100644 --- a/pyxtal/util.py +++ b/pyxtal/util.py @@ -381,7 +381,7 @@ def get_matching_coord(coord, ops, atol=1e-4): for i, d in enumerate(p._cif.data.values()): ops = p.get_symops(d) coord_to_species = OrderedDict() - d0 = {"_atom_site_label": [], + d0 = {"_atom_site_label": [], "_atom_site_fract_x": [], "_atom_site_fract_y": [], "_atom_site_fract_z": [], @@ -391,12 +391,12 @@ def get_matching_coord(coord, ops, atol=1e-4): symbol = p._parse_symbol(d["_atom_site_type_symbol"][i]) except KeyError: symbol = p._parse_symbol(d["_atom_site_label"][i]) - + el = get_el_sp(symbol) x = str2float(d["_atom_site_fract_x"][i]) y = str2float(d["_atom_site_fract_y"][i]) z = str2float(d["_atom_site_fract_z"][i]) - + coord = (x, y, z) match = get_matching_coord(coord, ops) if not match: @@ -409,7 +409,7 @@ def get_matching_coord(coord, ops, atol=1e-4): d.data['_atom_site_fract_x'] = d0['_atom_site_fract_x'] d.data['_atom_site_fract_y'] = d0['_atom_site_fract_y'] d.data['_atom_site_fract_z'] = d0['_atom_site_fract_z'] - + s = p._get_structure(d, primitive=False, symmetrized=False) return s @@ -488,10 +488,10 @@ def sort_by_dimer(atoms, N_mols, id=10, tol=4.0): atoms.set_positions(pos1) return atoms -def generate_wp_lib(spg_list, composition, - num_wp=(None, None), - num_fu=(None, None), - num_dof=(None, None), +def generate_wp_lib(spg_list, composition, + num_wp=(None, None), + num_fu=(None, None), + num_dof=(None, None), N_max=1000): """ Generate wps according to the composition constraint (e.g., SiO2) @@ -523,13 +523,13 @@ def generate_wp_lib(spg_list, composition, for sg in spg_list: g = Group(sg) lat_dof = g.get_lattice_dof() - # determine the upper and lower limit + # determine the upper and lower limit if min_fu is None: min_fu = max([int(len(g[-1])/min(composition)), 1]) if max_fu is None: max_fu = max([int(len(g[0])/max(composition)), 1]) count = 0 for i in range(max_fu, min_fu-1, -1): letters, _, wp_ids = g.list_wyckoff_combinations( - composition*i, max_wp=max_wp, + composition*i, max_wp=max_wp, min_wp=min_wp, Nmax=100000) for j, wp in enumerate(wp_ids): wp_dofs = 0 @@ -545,8 +545,23 @@ def generate_wp_lib(spg_list, composition, count += 1 if count >= N_max: break - return wps - + return wps + + +def reset_lammps_cell(atoms0): + """ + set the cell into lammps format + """ + from ase.calculators.lammpslib import convert_cell + + atoms = atoms0.copy() + mat, coord_transform = convert_cell(atoms0.cell) + if coord_transform is not None: + pos = np.dot(atoms0.positions, coord_transform.T) + atoms.set_cell(mat.T) + atoms.set_positions(pos) + return atoms + if __name__ == "__main__": from argparse import ArgumentParser