Rosa bulo/pdb and mol2

CHANGELOG.md

@@ -15,6 +15,8 @@ This changelog is effective from the 2025 releases.
 * MultiJob now supports generic Job types for the self.children attribute
 
 ### Changed
+* `Molecule.readmol2` can now read non-integer bond orders
+* `Molecule.writepdb` now writes default values for res, resnum, fix, and occ columns
 * `view` function uses stdin mode for AMSview, reducing overhead for image creation
 * `vasp_output_to_ams` now reads the MD time step from `POTIM` in the OUTCAR and labels molecular-dynamics runs (`IBRION = 0`) as such, instead of using a fixed default time step
 

src/scm/plams/mol/molecule.py

@@ -3184,6 +3184,8 @@ def readmol2(self, f: IO, **other: Any) -> None:
                 spl = line.split()
                 if len(spl) < 4:
                     raise FileError(f"readmol2: Error in {f.name} line {i+1}: not enough values in line")
+                if "." in spl[3]:
+                    spl[3] = f"{int(float(spl[3]))}"
                 try:
                     atom1 = self.atoms[int(spl[1]) - 1]
                     atom2 = self.atoms[int(spl[2]) - 1]
@@ -3278,15 +3280,25 @@ def writepdb(self, f: IO, **other: Any) -> None:
         pdb = PDBHandler()
         for i, at in enumerate(self.atoms):
             pdbatom = PDBAtom()
+            pdbatom.name = f"{at.symbol.upper():<2}"
+            pdbatom.name = "%-2s" % (at.symbol.upper())
             pdbatom.coords = at.coords
             pdbatom.element = at.symbol.upper()
+            pdbatom.res = "LIG"
+            pdbatom.resnum = "0"
+            pdbatom.occ = "1.00"
+            pdbatom.fix = "0.00"
             if "pdb" in at.properties:
                 if "res" in at.properties.pdb:
                     pdbatom.res = at.properties.pdb.res
                 if "resnum" in at.properties.pdb:
                     pdbatom.resnum = at.properties.pdb.resnum
                 if "name" in at.properties.pdb:
                     pdbatom.name = at.properties.pdb.name
+                if "occ" in at.properties.pdb:
+                    pdbatom.occ = at.properties.pdb.occ
+                if "fix" in at.properties.pdb:
+                    pdbatom.fix = at.properties.pdb.fix
             pdb.add_atom(pdbatom)
         if len(self.lattice) > 0:
             pdb.set_lattice(self.lattice)

unit_tests/test_molecule.py

@@ -1,4 +1,5 @@
 import os
+from io import StringIO
 import pytest
 from abc import ABC, abstractmethod
 
@@ -1113,6 +1114,42 @@ def test_write_multiple_molecules_to_pdb(pdb_folder, tmp_path):
             np.testing.assert_allclose(at.coords, at_ref.coords, atol=1e-08)
 
 
+def test_simple_pdb_writing(xyz_folder, tmp_path):
+    """
+    Test that a reasonable PDB format is written for a simple benzene molecule
+
+    The file should contain for each atom a default resname and resnum,
+    as well as default fix and occ columns.
+    """
+    benzene = Molecule(os.path.join(xyz_folder, "benzene.xyz"))
+    f = StringIO()
+    benzene.writepdb(f)
+    f.seek(0)
+    for i in range(2):
+        line = f.readline()
+
+    refline = "ATOM      1   C  LIG     0       1.194  -0.689   0.000  1.00  0.00           C  \n"
+    assert line == refline
+
+
+def test_mol2_format(xyz_folder, tmp_path):
+    """
+    Test that a Mol2 file can be written and read for a molecule with double and aromatic bonds
+    """
+    mol = Molecule(os.path.join(xyz_folder, "reactant2.xyz"))
+    mol.guess_bonds()
+
+    f = StringIO()
+    mol.writemol2(f)
+    f.seek(0)
+
+    mol2 = Molecule()
+    mol2.readmol2(f)
+
+    assert len(mol) == len(mol2)
+    assert len(mol.bonds) == len(mol2.bonds)
+
+
 def test_read_multiple_molecules_from_coskf(coskf_folder):
     """
     Test for COSKF reading