added test for shape consistency of masses, momenta and positions

classical/systems.py

@@ -161,6 +161,8 @@ class ManyParticlePhysicalSystem:
         self.mass = np.array(mass)
         self.cell = np.array(cell)
         self.ff = ff
+        assert self.mass.shape == self.q.shape == self.p.shape
+        assert self.mass.shape == self.q0.shape == self.p0.shape
         self.kinetic_energy = sum(p**2/(2*m) for p, m in zip(self.p, self.mass))
         self.potential_energy = self.ff.get_potential_energy('cartesian')
         self.hamiltonian = (self.potential_energy, self.kinetic_energy)