求助：压力耦合方法的选择求助

呆宝可乐

模拟体系包含12000CO2+5[Emim][Tf2N]+10Ls-45+30EtOH，想研究表面活性剂的自组装，我的问题是采用P-R压力耦合方法blow up严重，说是体系不平衡，而采用Berensden并没有崩，翻看论坛记录sob老师几年前有提到说大体系压力耦合方法不一定非要采用P-R方法，“

Berensden压浴从实际来讲完全没问题，有些人对理论过于吹毛求疵了”，但又有人说Berensden方法并不能反应真实的情况，那我想问，针对我这一体系组成，压浴是否可以采用Berensden方法？是否会对结果造成很大的影响。

sobereva

用Berensden完全没问题
就算你想用PR，也应当等到体系已经跑了一阵子，已经充分平衡了再改成PR

tjuptz

sobereva 发表于 2020-5-17 03:17
用Berensden完全没问题
就算你想用PR，也应当等到体系已经跑了一阵子，已经充分平衡了再改成PR

老师好，我前段时间学习ASE中的NPT模块（NH热浴+PR压浴）时遇到个参数问题：
pfactor:
A constant in the barostat differential equation. If a characteristic barostat timescale of ptime is desired, set pfactor to ptime^2 * B (where B is the Bulk Modulus). Set to None to disable the barostat. Typical metallic bulk moduli are of the order of 100 GPa or 0.6 eV/?^3. Useful parameter values:
The ttime and pfactor are quite critical[4], too small values may cause instabilites and/or wrong fluctuations in T / p. Too large values cause an oscillation which is slow to die. Good values for the characteristic times seem to be 25 fs for ttime, and 75 fs for ptime (used to calculate pfactor), at least for bulk copper with 15000-200000 atoms. But this is not well tested, it is IMPORTANT to monitor the temperature and stress/pressure fluctuations.
我查了下体香bulk modulus和compressibility是倒数关系，比如水的压缩系数是4.5*10^-5 bar^-1，这也是GROMACS里压浴的常用参数。
我又看了其源代码中pfactor和ptime的转换关系如下，有点迷惑，到底该怎么设置呢？
def _calculateconstants(self):
        "(Re)calculate some constants when pfactor, ttime or temperature have been changed."
        n = self._getnatoms()
        if self.ttime is None:
            self.tfact = 0.0
        else:
            self.tfact = 2.0 / (3 * n * self.temperature *self.ttime * self.ttime)
        if self.pfactor_given is None:
            self.pfact = 0.0
        else:
            self.pfact = 1.0 / (self.pfactor_given * linalg.det(self._getbox()))
            # self.pfact = 1.0/(n * self.temperature * self.ptime * self.ptime)
        self.desiredEkin = 1.5 * (n - 1) * self.temperature