Параметры PVT задавал по разному, последний вариант: термостат: с коррекцией энергии на краях распределения, давление: сила прикладывается к граничным молекулам, выделенным алгоритмом Convex Hull.
А термостат какой? И какие силы прикладываются, в соответствии с каким алгоритмом?
Кажется, разница у нас только в знаке перед вириалом.
Нет, там нет разницы со знаками, у меня минус в

ещё запихан. Собственно, я возмущался первой частью (

): Вы туда ставите заданную температуру, или моментальную, высчитанную из кинетической энергии? Записано, как первый вариант, что неверно, но если имелся в виду второй, то всё норм.
Объем считается динамически, по алгоритму Convex Hull (картинка).
Вообще без периодических условий, строго говоря, объём не определён, поэтому можно позволить себе пляски с бубнами. Например, Вы можете отмасштабировать этот объём каким-нть коэффициентом

, чтобы учесть, что молекулы воды -- это не точки, а размазанные в пространстве электронные плотности (по-сути добавить аналог solvent accessibility). Тогда флуктуации давления упадут по абсолютной величине.
Полно статей, где народ рапортует о том как хорошо совпадает диаграмма состояния для потенциалов TIP3P, TIP4P. Как же они считают давление?
Во-первых, Вы можете попробовать нормально определить объём, задав граничные условия. Во-вторых, Вас интересуют не моментальные величины, которые из-за флуктуаций могут быть какими угодно, а средние значения ансамбля, усреднённые по достаточно большому уравновешенному куску траектории. В-третьих, Вы можете всё то же самое проделать с нормальными баростатами, типа Берендсоновского.