计算化学公社

标题: 求助，为什么用NVT模拟气液界面上水蒸发时得到压力为负数 [打印本页]

作者
Author: BQHY 时间: 2019-9-28 19:08
标题: 求助，为什么用NVT模拟气液界面上水蒸发时得到压力为负数
从官网下的in文件，得到的表面张力也不对。恳求各位看下哪出了问题。。

units real
atom_style full
read_data data.singleTIP4P-2005
include forcefield.TIP4P-2005
pair_modify tail no

replicate 20 20 20

variable xLo equal -0.5*lx
variable xHi equal 0.5*lx
variable yLo equal -0.5*ly
variable yHi equal 0.5*ly
variable zLo equal -0.5*lz
variable zHi equal 0.5*lz

# Recenter system coords about the origin (0,0,0):
displace_atoms all move ${xLo} ${yLo} ${zLo} units box

# Adjust box dimensions:
change_box all x final ${xLo} ${xHi} units box
change_box all y final ${yLo} ${yHi} units box
change_box all z final ${zLo} ${zHi} units box

# Increase z edge to yield "Vacuum|Water|Vacuum" system:
variable zLoNew equal -2.5*lz
variable zHiNew equal 2.5*lz
change_box all z final ${zLoNew} ${zHiNew} units box

variable Nblock equal 10000
variable Nrun equal 2.0*${Nblock}
# variable Ndump equal ${Nrun}/2
variable Ne equal 100
variable Nr equal ${Nblock}/${Ne}
variable Nf equal ${Nrun}
variable Dz equal 10

variable A_in_m equal 1e-10 # Angstrom in meter
variable atm_in_Pa equal 101325 # note: 1 Pa = 1 N/m^2
variable N_in_mN equal 1e3 # Newton in milliNewton
variable convFac equal ${A_in_m}*${atm_in_Pa}*${N_in_mN}

variable Text equal 340

group hydrogen type 1
group oxygen type 2

velocity all create ${Text} 12345

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

timestep 2.0

# fix loadBalance all balance 20 z 5 1.05

fix constrain all shake 1.0e-4 100 0 b 1 a 1
fix integrate all nvt temp ${Text} ${Text} 200.0
fix removeMomentum all momentum 1 linear 1 1 1

compute T all temp
fix TAve all ave/time ${Ne} ${Nr} ${Nf} c_T

variable P equal press
fix PAve all ave/time ${Ne} ${Nr} ${Nf} v_P

variable xPress equal c_thermo_press[1]
variable yPress equal c_thermo_press[2]
variable zPress equal c_thermo_press[3]

# Evaluate and average surface tension in mN/m:
variable st equal 0.5*lz*(v_zPress-0.5*(v_xPress+v_yPress))*${convFac}
fix st all ave/time ${Ne} ${Nr} ${Nf} v_st

variable xyArea equal lx*ly

thermo_style custom step temp f_TAve press f_PAve v_xyArea lz f_st
thermo_modify flush yes
thermo ${Nrun}

compute A all chunk/atom bin/1d z center ${Dz} units box
fix massDensA all ave/chunk ${Ne} ${Nr} ${Nf} A density/mass &
file massDensA.zProfile

dump trj all atom 100 wat.lammpstrj
dump_modify trj scale no sort id
dump 123 all custom 200 dump.h2o.lv.* id type x y z vx vy vz c_A

run ${Nrun}

作者
Author: 灰天皮匠 时间: 2019-9-29 10:13
本帖最后由灰天皮匠于 2019-9-29 10:17 编辑

这个in文件看不太懂，倒是可以提一些建议。
1、看一下力场参数是不是不合适
2、搞清楚里面命令含义以及如何设置。就像timestep=2.0和nvt控温dt=200，别人是不是这样模拟的
3、输出文件dump中加入总能量看看
图片是lammps手册中关于压强公式的描述

作者
Author: BQHY 时间: 2019-9-29 19:38

灰天皮匠发表于 2019-9-29 10:13
这个in文件看不太懂，倒是可以提一些建议。
1、看一下力场参数是不是不合适
2、搞清楚里面命令含义以及如 ...

好的，我再看看，非常感谢。

作者
Author: Picardo 时间: 2019-10-1 09:40
这个问题我没有做过，除了楼上上的建议我还有一个额外的看法，楼主有没有用其他系综看看计算结果，相变过程达到平衡态，两相的压力要相等，为何不用NPT系综呢

作者
Author: 961106 时间: 2019-12-17 20:38

BQHY 发表于 2019-9-29 19:38
好的，我再看看，非常感谢。

老师您好，请问可以提供一下tip4p/2005.itp的水模型文件嘛？我最近也想用这个水模型，但是不知道从哪里可以找，谢谢老师！

作者
Author: BQHY 时间: 2019-12-21 16:19

961106 发表于 2019-12-17 20:38
老师您好，请问可以提供一下tip4p/2005.itp的水模型文件嘛？我最近也想用这个水模型，但是不知道从哪里可 ...

我是在这个网站下的，你可以看看。http://www.orsi.sems.qmul.ac.uk/downloads.html

作者
Author: naoki 时间: 2019-12-21 16:38

Picardo 发表于 2019-10-1 09:40
这个问题我没有做过，除了楼上上的建议我还有一个额外的看法，楼主有没有用其他系综看看计算结果，相变过程 ...

NPT会不会直接把真空区压没了

作者
Author: 961106 时间: 2020-1-3 09:36

BQHY 发表于 2019-12-21 16:19
我是在这个网站下的，你可以看看。http://www.orsi.sems.qmul.ac.uk/downloads.html

谢谢

欢迎光临计算化学公社 (http://bbs.keinsci.com/)