Add a test case for parsing a DRUDE PSF and fix some issues with it

parmed/charmm/psf.py

@@ -156,9 +156,9 @@ def _parse_psf_section(cls, psf):
                 title, pointers = cls._parse_psf_title_line(line)
                 line = psf.readline().strip()
         data = []
-        if title == 'NATOM' or title == 'NTITLE':
-            # Store these two sections as strings (ATOM section we will parse
-            # later). The rest of the sections are integer pointers
+        if title in {'NATOM', 'NTITLE', 'NUMLP NUMLPH', 'NUMANISO'}:
+            # Store these two sections as strings to be parsed later.
+            # The rest of the sections are integer pointers
             while line:
                 data.append(line)
                 line = psf.readline().strip()
@@ -202,7 +202,7 @@ def __init__(self, psf_name=None):
             psfsections = _ZeroDict()
             while True:
                 try:
-                    sec, ptr, data = CharmmPsfFile._parse_psf_section(fileobj)
+                    sec, ptr, data = self._parse_psf_section(fileobj)
                 except _FileEOF:
                     break
                 psfsections[sec] = (ptr, data)
@@ -244,7 +244,7 @@ def __init__(self, psf_name=None):
                 if is_drude and i > 1 and drude_alpha_thole[-2] != (0, 0):
                     my_alpha, my_thole = drude_alpha_thole[-2]
                     atom = DrudeAtom(name=name, type=attype, charge=charge, mass=mass,
-                                     atomic_number=0, alpha=my_alpha, thole=my_thole, drude_type=atom_type)
+                                     atomic_number=0, alpha=my_alpha, thole=my_thole, drude_type=attype)
                 elif name.startswith('LP') and mass == 0:
                     atom = ExtraPoint(name=name, type=attype, charge=charge, mass=mass, atomic_number=0)
                 else:
@@ -400,7 +400,7 @@ def small_to_zero(num: float) -> float:
             return num
 
         ang *= DEG_TO_RAD
-        dihed *= DEG_TO_RAD
+        dihed = (180 - dihed) * DEG_TO_RAD
         if dist > 0:
             x_weights = [-1.0, 0.0, 1.0]
         elif dist < 0:

test/test_parmed_charmm.py

@@ -11,7 +11,7 @@
 from parmed.charmm import charmmcrds, parameters, psf
 from parmed.charmm._charmmfile import CharmmFile, CharmmStreamFile
 from parmed import exceptions, topologyobjects as to, load_file, ParameterSet
-from parmed.topologyobjects import BondType, AngleType, DihedralType, DihedralTypeList
+from parmed.topologyobjects import BondType, AngleType, DihedralType, DihedralTypeList, DrudeAtom, ExtraPoint
 import parmed.unit as u
 import random
 import unittest
@@ -1181,3 +1181,49 @@ def test_write_xplor(self):
         parm.save(fn, overwrite=True)
         cpsf = pmd.load_file(fn)
         self.assertIn('XPLOR', cpsf.flags)
+
+class TestCharmmDrudeSystems(unittest.TestCase):
+    """ Tests processing of DRUDE systems """
+
+    def test_drude_parsing_na(self):
+        """ Test parsing a trinucleotide with Drude particles and lone pairs """
+        system = psf.CharmmPsfFile(get_fn("cyt-gua-cyt.psf"))
+        self.assertEqual(len(system.atoms), 176)
+
+        # Check that we have the right number of DRUDE particles
+        self.assertEqual(sum(isinstance(atom, DrudeAtom) for atom in system.atoms), 59)
+        self.assertTrue(all(atom.type == "DRUD" for atom in system.atoms if isinstance(atom, DrudeAtom)))
+
+        # Check that we have the right number of lone pairs
+        self.assertEqual(sum(isinstance(atom, ExtraPoint) for atom in system.atoms), 23)
+
+        # Check that the parameters match expectation from a known implementation
+        ep1 = system.atoms[7]
+        origin_weights, x_weights, y_weights, local_position = ep1.frame.get_weights()
+        self.assertAlmostEqual(local_position[0], -0.11970705016398403)
+        self.assertAlmostEqual(local_position[1], 0.005739964140425085)
+        self.assertAlmostEqual(local_position[2], 0.32884232536689074)
+        self.assertEqual(x_weights, [-1, 0, 1])
+        self.assertEqual(y_weights, [0, -1, 1])
+        self.assertEqual(origin_weights, [1, 0, 0])
+
+        # Now try a virtual site with a negative distance
+        ep2 = system.atoms[22]
+        origin_weights, x_weights, y_weights, local_position = ep2.frame.get_weights()
+        self.assertAlmostEqual(local_position[0], -0.34999999466919514)
+        self.assertAlmostEqual(local_position[1], 6.108652257915055e-05)
+        self.assertAlmostEqual(local_position[2], 1.0661609584247726e-08)
+        self.assertEqual(x_weights, [-1, 0.5, 0.5])
+        self.assertEqual(y_weights, [0, -1, 1])
+        self.assertEqual(origin_weights, [1, 0, 0])
+
+        # Check some Drude particle parameters
+        drude1 = system.atoms[2]
+        self.assertIsInstance(drude1, DrudeAtom)
+        self.assertEqual(drude1.alpha, -1.028)
+        self.assertEqual(drude1.thole, 1.3)
+
+        drude2 = system.atoms[-1]
+        self.assertIsInstance(drude2, DrudeAtom)
+        self.assertEqual(drude2.alpha, -0.157)
+        self.assertEqual(drude2.thole, 1.3)