Bug fix of bond_order_system dump to sdf file

README.md

@@ -1,4 +1,4 @@
-**dpdata** is a python package for manipulating DeePMD-kit, VASP, LAMMPS data formats.
+**dpdata** is a python package for manipulating data formats of software in computational science, including DeePMD-kit, VASP, LAMMPS, GROMACS, Gaussian.
 dpdata only works with python 3.x.
 
 
@@ -116,7 +116,7 @@ xyz_multi_systems.to_deepmd_raw('./my_deepmd_data/')
 ```
 
 You may also use the following code to parse muti-system:
-```
+```python
 from dpdata import LabeledSystem,MultiSystems
 from glob import glob
 """

dpdata/bond_order_system.py

@@ -4,6 +4,7 @@
 import dpdata.rdkit.utils
 from dpdata.rdkit.sanitize import Sanitizer, SanitizeError
 from copy import deepcopy
+from rdkit.Chem import Conformer
 # import dpdata.rdkit.mol2
 
 def check_BondOrderSystem(data):
@@ -95,6 +96,12 @@ def from_fmt_obj(self, fmtobj, file_name, **kwargs):
         return self
 
     def to_fmt_obj(self, fmtobj, *args, **kwargs):
+        self.rdkit_mol.RemoveAllConformers()
+        for ii in range(self.get_nframes()):
+            conf = Conformer()
+            for idx in range(self.get_natoms()):
+                conf.SetAtomPosition(idx, self.data["coords"][ii][idx])
+            self.rdkit_mol.AddConformer(conf, assignId=True)
         return fmtobj.to_bond_order_system(self.data, self.rdkit_mol, *args, **kwargs)
 
     def __repr__(self):
@@ -151,7 +158,8 @@ def copy(self):
         self.__class__(data=deepcopy(self.data),
                        rdkit_mol=new_mol)
 
-    # def __add__(self, other):
+    def __add__(self, other):
+        raise NotImplementedError("magic method '+' has not been implemented on BondOrderSystem")
     #     '''
     #         magic method "+" operation
     #     '''

dpdata/plugins/rdkit.py

@@ -13,8 +13,8 @@ def from_bond_order_system(self, file_name, **kwargs):
         return rdkit.Chem.MolFromMolFile(file_name, sanitize=False, removeHs=False)
 
 
-    def to_bond_order_system(self, mol, data, file_name, frame_idx=0, **kwargs):
-        assert (frame_idx < self.get_nframes())
+    def to_bond_order_system(self, data, mol, file_name, frame_idx=0, **kwargs):
+        assert (frame_idx < mol.GetNumConformers())
         rdkit.Chem.MolToMolFile(mol, file_name, confId=frame_idx)
 
 
@@ -32,12 +32,12 @@ def from_bond_order_system(self, file_name, **kwargs):
             mol = mols[0]
         return mol
 
-    def to_bond_order_system(self, mol, data, file_name, frame_idx=-1, **kwargs):
+    def to_bond_order_system(self, data, mol, file_name, frame_idx=-1, **kwargs):
         sdf_writer = rdkit.Chem.SDWriter(file_name)
         if frame_idx == -1:
-            for ii in range(data['coords'].shape[0]):
+            for ii in range(mol.GetNumConformers()):
                 sdf_writer.write(mol, confId=ii)
         else:
-            assert (frame_idx < data['coords'].shape[0])
+            assert (frame_idx < mol.GetNumConformers())
             sdf_writer.write(mol, confId=frame_idx)
         sdf_writer.close()

dpdata/rdkit/utils.py

@@ -63,7 +63,7 @@ def system_data_to_mol(data):
     mol = mol_ed.GetMol()
     # set formal charges
     for idx, atom in enumerate(mol.GetAtoms()):
-        atom.SetFormalCharge(data['formal_charges'][idx])
+        atom.SetFormalCharge(int(data['formal_charges'][idx]))
     # set mol name
     if '_name' in list(data.keys()):
         mol.SetProp("_Name", data['_name'])

tests/test_bond_order_system.py

@@ -6,6 +6,7 @@
 from rdkit.Chem import AllChem
 import shutil
 import numpy as np
+from copy import deepcopy
 
 class TestBondOrderSystem(unittest.TestCase):
 
@@ -78,6 +79,16 @@ def test_dump_to_deepmd_npy(self):
                 self.assertEqual(syst['bonds'][bond_idx][ii], bonds[bond_idx][ii])
         shutil.rmtree("bond_order/methane")
 
+    def test_dump_to_sdf_file(self):
+        s1 = dpdata.BondOrderSystem("bond_order/methane.sdf", fmt="sdf")
+        s2 = deepcopy(s1)
+        s2.data["coords"] += 1.0
+        s2.to_sdf_file("bond_order/test.sdf")
+
+        nsyst = dpdata.BondOrderSystem("bond_order/test.sdf", fmt="sdf")
+        self.assertEqual(nsyst["coords"][0, 0, 0] - s1["coords"][0, 0, 0], 1.0)
+        os.remove("bond_order/test.sdf")
+
     def test_sanitize_mol_obabel(self):
         cnt = 0
         for sdf_file in glob.glob("bond_order/refined-set-ligands/obabel/*sdf"):