チュートリアル11:界面張力計算

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

サンプリング

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

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

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

なお、新たに気液界面の系を作成する場合は、平衡化時に壁を置くなどinputファイルを変更する必要がある。（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を参照）

また、圧力テンソルの計算においては拘束力計算の収束判定条件をデフォルトの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

と実行すると、以下のようなプロフィールを得る。

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:界面張力計算

サンプリング

サンプリング結果

案内メニュー

検索