MDAnalysis / mdanalysis

MDAnalysis is a Python library to analyze molecular dynamics simulations.
https://mdanalysis.org
Other
1.33k stars 653 forks source link

MDAnalysis.analysis.dielectric fails when dealing with tip4p #4119

Open VOD555 opened 1 year ago

VOD555 commented 1 year ago

Expected behavior

MDAnalysis.analysis.dielectric should give the dielectric of the system.

Actual behavior

  1. Given make_whole=True, it can't handle the 0 mass virtual sites as discuessed in #3168
  2. Given make_whole=False, it still can't work well, as the dummy site are not included in the same fragment as the O and H atoms as there's no bond defiend between the dummy site and other atoms in the water model, resulting in non-neutral total charge of the fragments. And mdanalysis just thinks they are free charges. This is related to #1954

My current way to solve this is to make a fake .tpr file that contains a bond between the virtual site and the oxygen.

Code to reproduce the behavior

import MDAnalysis as mda
from MDAnalysis.analysis import dielectric 

u = mda.Universe('md.tpr', 'md.xtc')

diel = dielectric.DielectricConstant(u.atoms, temperature=298, make_whole=False, start=5000)
diel.run()

Current version of MDAnalysis

orbeckst commented 1 year ago

@VOD555 , MAICoS has some tricks related to TIP4P https://maicos.readthedocs.io/en/main/examples/dielectric-profiles.html#tip4p-water-and-molecules-with-virtual-sites — @PicoCentauri can probably tell you more.