add a function to generate forcefield from ost, experimental

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")

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()

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

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