Hirano: ページの作成:「界面張力の計算を行う時、平衡化については他の方法と同様である。従ってここでは、サンプリング(インプットとしてどの変…」

2026-05-26T02:37:48Z

ページの作成:「界面張力の計算を行う時、平衡化については他の方法と同様である。従ってここでは、サンプリング(インプットとしてどの変…」

新規ページ

界面張力の計算を行う時、平衡化については他の方法と同様である。従ってここでは、サンプリング(インプットとしてどの変数を指定するか)と、それによって得られたアウトプットファイルへの処理について述べることとする。

== サンプリング ==
ここでは、平衡化によって80トラジェクトリ―のファイルが得られたとき、そこからどのように界面張力を計算していくのかを水の気液界面を例にとって述べることとする。
また、この計算におけるセルサイズなどは以下のようになっている。

[[ファイル:Water_set.png|500px]]

計算を行う際は、tutorial11内部にあるファイルを一旦FreeFlex外部にコピーしてから行うこと。

今回の計算では、80個の系を計算し、それぞれを平均して圧力を求める。以下にサブミットファイルの例を示す。

なお、新たに気液界面の系を作成する場合は、平衡化時に壁を置くなどinputファイルを変更する必要がある。（[[チュートリアル10:溶液界面のグランドカノニカルMD|Tutorial10]]を参照）
=== sub_fep.sh ================================================================================
#!/bin/bash
#PBS -N tutorial_11
#PBS -l select=1:ncpus=16:mpiprocs=16:ompthreads=1
#PBS -l walltime=80:00:00
#PBS -j oe

###### ユーザー指定の変数の default 以外を設定する。###########################
RNODE=****00
RDIR=
RFILE=
WDIR=/work/(USER)/$$
SAVE=

echo -n 'START: '; date
###### initialization #####################################################
###### [1] mpd ######
touch /tmp/node_$$
trap 'rm /tmp/node_$$' 1 2 3 15
uniq $PBS_NODEFILE > /tmp/node_$$
NPROCS=`cat $PBS_NODEFILE | wc -l`
NHOSTS=`cat /tmp/node_$$ | wc -l`
for i in `cat /tmp/node_$$`
do
rsh $i mpdcleanup >/dev/null
done
mpdboot -f /tmp/node_$$ -n $NHOSTS -r rsh

###### [2] files ######
echo 'RNODE: '${RNODE:='fep'} # default は、'fep'になる。
echo 'RDIR: '${RDIR:=$PBS_O_WORKDIR} # defalut は、サブミットした場所。
echo 'RFILE: '${RFILE:="*"} # default は、すべて。
echo 'WDIR: '${WDIR:=$HOME/$$} # default は、$HOME/$$
echo 'SAVE: '${SAVE:=0} # default は、0
echo 'PBS_NODE: '`cat /tmp/node_$$`

test -e $WDIR || mkdir -p $WDIR # WDIR をつくる。
cd $WDIR
touch touch_$$
if [ "`rsh $RNODE ls $RDIR/touch_$$ 2>/dev/null`" == "$RDIR/touch_$$" ]; then
IDENTICAL=0 # RDIR と WDIR が同一のとき、何もしない。
else
IDENTICAL=1 # RDIR と WDIR が別ディレクトリのとき、RNODE の入力ファイルを WDIR にコピー。
rsh $RNODE "cd $RDIR; tar czf - $RFILE" | tar xzf -
fi
rm touch_$$

####### execution ###########################################################
mpiexec -machinefile $PBS_NODEFILE -n $NPROCS /home/(USER)/FreeFlex/FreeFlex.exe input_sam.nml
###############################################################################
...
===============================================================================================

こういったサブミットファイルを作業ディレクトリ上で
> qsub sub_fep.sh
とすれば、計算が行われる。

この計算は計算機にもよるが大体３日ほどかかるので、必要に応じてompthreadを変更するとよい。

サンプリングの際のnmlファイル(input_sam.nml)を以下に示す。
=== input_sam.nml =============================================================================
&FreeFlex
（略）
calculate_press = .true.
/

&MONITOR
interval = 1
file="press_sam.output"
condition = "PRESS"
overwrite = .false.
/

&MONITOR
interval = 1
file="pressCON_sam.output"
condition = "PRESSCON"
overwrite = .false.
/

&MONITOR
interval = 1
file="press_z_sam.output"
condition = "PRESS_z"
overwrite = .true.
/

&MONITOR
interval = 1
file="pressCON_z_sam.output"
condition = "PRESSCON_z"
overwrite = .true.
/

&constraint
tolerance = 1.0d-18
/
===============================================================================================
ここで、FreeFlexのネームリストとして、calculate_press = .true.とすることで圧力テンソルの計算を行う。

圧力テンソルの計算結果を出力する変数としては、PRESSS、PRESSCON、PRESS_ｚ、PRESSCON_ｚがあり、これらを指定する必要がある。この内PRESS_z、PRESSCON_zについては、overwrite=.true.としないと毎ステップ各ビンの圧力テンソルが出力されるので、これを.true.とすることでステップについて平均の値を出力することを推奨する。またこれらの変数の説明についてはネームリストのページを見てほしい。([[Monitor System|&MONITOR]]を参照）

また、圧力テンソルの計算においては拘束力計算の収束判定条件をデフォルトの1.0d-6 angstromsではなく、1.0d-8 angstromにしなければならない。

== サンプリング結果 ==
MONITOR変数の指定によって界面に対して垂直な方向に沿った圧力テンソルの出力として、以下のようなファイルが得られる。
========press_z_sam.output========================================================================================================================================================
z, P_K[xx,yy,zz], P_U [xx, yy, zz, xy, yz, zx]
-0.4950000000E-08 0.2805166548E+05 0.4418980353E+05 0.5998336912E+05 -0.6615222747E+04 -0.9190008810E+04 0.7859812530E+04 0.1358711871E+04 -0.9372980979E+03 -0.1105971529E+04
-0.4850000000E-08 0.3284841367E+05 0.4179705637E+05 0.5919083945E+05 0.3812638440E+04 0.5554852780E+04 -0.1714722034E+04 -0.1484721547E+04 0.1834053745E+04 0.8854751769E+03
・・・・・・・・・・・・・・・・・(略)
==================================================================================================================================================================================

========pressCON_z_sam.output==============================================================================================================
z, P_CON [xx, yy, zz, xy, yz, zx]
　-0.4950000000E-08 -0.2283145371E+05 -0.5107960470E+05 -0.4567547172E+05 -0.1930124471E+05 0.9124466958E+04 -0.1273752864E+05
-0.4850000000E-08 -0.2719560202E+05 -0.4385404221E+05 -0.4487171913E+05 -0.1161909582E+05 0.1979643347E+04 -0.1472059856E+05
・・・・・・・・・・・・・・・・・(略)
===========================================================================================================================================
トラジェクトリを80個並列で流した場合、このようなファイルも80個分出力される。これら80個のデータを平均化するには、以下のようなファイルを作り実行すればよい。
========ave_z.f90==============================================================================================================
program main
integer,parameter :: nfile=80
integer,parameter :: nline=100
character(len=4) :: cfile
real*8 :: z(1:nline), p(1:15,1:nline)
real*8 :: p_data(1:6,1:nline,1:nfile), p_ave(1:6,1:nline), p_diff(1:6,1:nline,1:nfile), error(1:6,1:nline)
real*8 :: z_, p_(1:15)
integer :: ifile,i

z = 0.d0
p = 0.d0
do ifile=0,nfile-1
write(cfile,'(i4.4)') ifile
write(*,*) cfile

open(90,file='../'//cfile//'/press_z_sam.output')
read(90,*)
do iline=1,nline
read(90,*) z_, p_(1:9)
z(iline) = z(iline) + z_
p(1:9,iline) = p(1:9,iline) + p_(1:9)
p_data(1:3,iline,ifile) = p_data(1:3,iline,ifile) + p_(1:3) +p_(4:6)
p_data(4:6,iline,ifile) = p_data(4:6,iline,ifile) + p_(7:9)
end do
close(90)

open(80,file='../'//cfile//'/pressCON_z_sam.output')
read(80,*)
do iline=1,nline
read(80,*) z_, p_(10:15)
p(10:15,iline) = p(10:15,iline) + p_(10:15)
p_data(1:6,iline,ifile) = p_data(1:6,iline,ifile) + p_(10:15)
end do
close(80)

end do
z = z/dble(nfile)
p = p/dble(nfile)

!!!error計算
p_ave(1:3,1:nline) = p(1:3,1:nline) + p(4:6,1:nline) +p(10:12,1:nline)
p_ave(4:6,1:nline) = p(7:9,1:nline) + p(13:15,1:nline)
do i = 1,nfile
p_diff(1:6,1:nline,i) = p_data(1:6,1:nline,i) - p_ave(1:6,1:nline)
end do
error(1:6,1:nline) = sqrt(sum(p_diff(1:6,1:nline,1:nfile)**2,3)/dble(nfile))

open(100,file='press_z_sam.dat')
do iline= 1,nline
===============================================================================================================================
これによって得られたpress_z_sam.datを以下のpz.pltとgnuplotを用いてプロットすれば良い。
========pz.plt=====================================================================================
reset

set encoding iso
set terminal png font "Times"

set out 'pol3_pressz.png'

set key font ",13"
set tics font ",13"
set xlabel font ",17"
set ylabel font ",17"
set grid

set xlabel "{/:Italic z} (\305)" offset 0.0,0.2

set ylabel "{/:Italic P} (MPa)" offset 0.0,0.0
set key top outside

a=10e-6
angs=1.e+10
dLz=9.91*10e-11

#P_K(運動エネルギー由来の圧力テンソル)
#plot 'press_z_sam.dat' u ($1*angs):($2*a) w l lw 2 lc rgb "dark-orange" title "xx",\
# 'press_z_sam.dat' u ($1*angs):($3*a) w l lw 2 lc rgb "orange" title "yy",\
# 'press_z_sam.dat' u ($1*angs):($4*a) w l lw 2 lc rgb "red" title "zz",\

#P_U(相互作用由来の圧力テンソル)
#plot 'press_z_sam.dat' u ($1*angs):($5*a) w l lw 2 lc rgb "dark-orange" title "xx",\
# 'press_z_sam.dat' u ($1*angs):($6*a) w l lw 2 lc rgb "orange" title "yy",\
# 'press_z_sam.dat' u ($1*angs):($7*a) w l lw 2 lc rgb "red" title "zz",\
# 'press_z_sam.dat' u ($1*angs):($8*a) w l lw 2 lc rgb "dark-blue" title "xy",\
# 'press_z_sam.dat' u ($1*angs):($9*a) w l lw 2 lc rgb "blue" title "yz",\
# 'press_z_sam.dat' u ($1*angs):($10*a) w l lw 2 lc rgb "skyblue" title "zx"

#P_CON(拘束力由来の圧力テンソル)
#plot 'press_z_sam.dat' u ($1*angs):($11*a) w l lw 2 lc rgb "dark-orange" title "xx",\
# 'press_z_sam.dat' u ($1*angs):($12*a) w l lw 2 lc rgb "orange" title "yy",\
# 'press_z_sam.dat' u ($1*angs):($13*a) w l lw 2 lc rgb "red" title "zz",\
# 'press_z_sam.dat' u ($1*angs):($14*a) w l lw 2 lc rgb "dark-blue" title "xy",\
# 'press_z_sam.dat' u ($1*angs):($15*a) w l lw 2 lc rgb "blue" title "yz",\
# 'press_z_sam.dat' u ($1*angs):($16*a) w l lw 2 lc rgb "skyblue" title "zx"

#P(圧力テンソル)
plot 'press_z_sam.dat' u ($1*angs):($2*a+$5*a+$11*a) w l lw 2 lc rgb "dark-orange" title "xx",\
'press_z_sam.dat' u ($1*angs):($3*a+$6*a+$12*a) w l lw 2 lc rgb "orange" title "yy",\
'press_z_sam.dat' u ($1*angs):($4*a+$7*a+$13*a) w l lw 2 lc rgb "red" title "zz",\
'press_z_sam.dat' u ($1*angs):($8*a+$14*a) w l lw 2 lc rgb "dark-blue" title "xy",\
'press_z_sam.dat' u ($1*angs):($9*a+$15*a) w l lw 2 lc rgb "blue" title "yz",\
'press_z_sam.dat' u ($1*angs):($10*a+$16*a) w l lw 2 lc rgb "skyblue" title "zx"

##errorの表示
#plot 'press_z_sam.dat' u ($1*angs):($2*a+$5*a+$11*a):($17*a) with errorbars lc rgb "dark-orange" title
"xx",\
# 'press_z_sam.dat' u ($1*angs):($3*a+$6*a+$12*a):($18*a) with errorbars lc rgb "orange" title
"yy",\
# 'press_z_sam.dat' u ($1*angs):($4*a+$7*a+$13*a):($19*a) with errorbars lc rgb "red" title "zz",\
# 'press_z_sam.dat' u ($1*angs):($8*a+$14*a):($20*a) with errorbars lc rgb "dark-blue" title "xy",\
# 'press_z_sam.dat' u ($1*angs):($9*a+$15*a):($21*a) with errorbars lc rgb "blue" title "yz",\
# 'press_z_sam.dat' u ($1*angs):($10*a+$16*a):($22*a) with errorbars lc rgb "skyblue" title "zx"
===================================================================================================

これを
>gnuplot pz.plt
と実行すると、以下のようなプロフィールを得る。

[[ファイル:pressz.png|500px]]

ここで、プロットした結果において非対角成分(xy,xz,xz)が0に近くならない場合、何かしら異常があると思われる。

（例として、平衡化、サンプリング時間不足、初期構造が悪いなど）

また、この圧力テンソルのプロフィールについて

法線方向と接線方向の差の積分値を計算する以下のプログラムを実行すれば界面張力を得る。
===ave_g.f90===================================================================================================================
program main
implicit none
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer,parameter :: nfile=80
integer,parameter :: nbin=100
integer,parameter :: skip = 1
character(len=4) :: cfile
real*8 :: z, P_K(1:3), P_U(1:6), P_CON(1:6), P(1:6)
real*8 :: G, total_G(1:nfile), total_G2(1:16)
real*8 :: gamma, gamma2, Var_gamma
integer :: ibin,ifile,i,ifile_
real*8,parameter :: Lz = 100.0d-10 !!!z方向のセル長
real*8,parameter :: dz = 100.0d-12 !!!一つのビンの長さ(=Lz/nbin)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
total_G =0.0d0
total_G2 = 0.0d0
do ifile=0,nfile-1
write(cfile,'(i4.4)') ifile
write(*,*) cfile
ifile_ = ifile + 1

P = 0.0d0
G = 0.0d0
open(90,file='../'//cfile//'/press_z_sam.output',status='old')
!!!読み飛ばす
do i = 1,skip
read(90,*)
end do

open(80,file='../'//cfile//'/pressCON_z_sam.output',status='old')
!!!読み飛ばす
do i = 1,skip
read(80,*)
end do

do ibin=1,nbin
read(90,*) z,P_K(1:3),P_U(1:6)
read(80,*) z,P_CON(1:6)
P(1:3) = P_K(1:3)+P_U(1:3)+P_CON(1:3)
P(4:6) = P_U(4:6)+P_CON(4:6)
G = (P(3)-0.5d0*(P(1)+P(2)))*dz/2.0d0 !１つのセルにつき二つの界面で計算しているので2で割る。
total_G(ifile_) = total_G(ifile_) + G
end do
close(90)
close(80)

end do

write(*,*) 'total_G'
write(*,'(10e15.7)') total_G(1:nfile)

gamma = sum(total_G)/dble(nfile)
total_G(1:nfile) = total_G(1:nfile) - gamma
Var_gamma = sum(total_G(1:nfile)**2)/dble(nfile)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

write(*,*) 'surface tension (N/m)'
write(*,'(10e15.7)') gamma
write(*,*) 'deviation'
write(*,'(10e15.7)') sqrt(Var_gamma)

end program main
===============================================================================================================================

これを実行した結果は以下のようになった。
===================================================================================
total_G
0.4822124E-01 0.5462081E-01 0.5266729E-01 0.5506123E-01 0.5436713E-01 0.4681410E-01 0.5064201E-01
0.4967185E-01 0.4962646E-01 0.5293858E-01
0.4761314E-01 0.5303573E-01 0.5199197E-01 0.5304727E-01 0.5058411E-01 0.4977264E-01 0.4680782E-01
0.5084134E-01 0.5217295E-01 0.5207262E-01
0.5490097E-01 0.5208478E-01 0.4912948E-01 0.5185132E-01 0.5437104E-01 0.5169138E-01 0.5084999E-01
0.4665278E-01 0.4949034E-01 0.5567725E-01
0.5376276E-01 0.5171501E-01 0.5185965E-01 0.5044360E-01 0.5070019E-01 0.5199373E-01 0.5094507E-01
0.4988094E-01 0.4810289E-01 0.5059433E-01
0.5363357E-01 0.4834083E-01 0.5496372E-01 0.5069184E-01 0.5198555E-01 0.5078426E-01 0.5087554E-01
0.4956766E-01 0.5099645E-01 0.5549398E-01
0.5000711E-01 0.5323785E-01 0.5171341E-01 0.5418007E-01 0.4855018E-01 0.5104953E-01 0.5220878E-01
0.5099397E-01 0.5179496E-01 0.5348767E-01
0.5072900E-01 0.5494645E-01 0.4968401E-01 0.5129792E-01 0.5289388E-01 0.5010409E-01 0.4777386E-01
0.5578830E-01 0.5246478E-01 0.5441176E-01
0.5195844E-01 0.5194620E-01 0.4918779E-01 0.5375268E-01 0.5264301E-01 0.4842055E-01 0.5304154E-01
0.5348141E-01 0.5046393E-01 0.5474954E-01
surface tension (N/m)
0.5154455E-01
derivation
0.2209831E-02

===================================================================================

ここでdeviationは80個のデータにおける標準偏差であり、標準誤差はこれの約√80分の一である。

チュートリアル11:界面張力計算 - 版の履歴

Hirano: ページの作成:「界面張力の計算を行う時、平衡化については他の方法と同様である。従ってここでは、サンプリング(インプットとしてどの変…」