튜토리얼 1: 가스 하이드레이트 구조 제작 및 시뮬레이션 해보기

::목표::

1. genice source code를 활용해 Clathrate hydrate의 구조들을 제작해볼 수 있습니다.
2. GROMACS 구동에 필요한 기본 파일 타입들을 학습합니다.
3. 생성한 분자 좌표 파일과 GROMACS 내장 Force fields를 활용해 간단히 MD 시뮬레이션을 진행해봅니다.
4. 진행된 시뮬레이션을 간단히 분석해볼 수 있습니다.

Table of Contents:

1. GenIce 소개

(1) Repository 간단 소개

Github를 둘러보면, 다양한 사람들이 작성해놓은 다양한 코드들이 존재하는 것을 알 수 있습니다.

그 중에는 하이드레이트 시뮬레이션에 특히 유용한 소스 코드들도 즐비하는데, 이 튜토리얼에서는 GenIce repository를 소개하지 않을 수가 없습니다.

Okayama 대학교의 Masakazu Matsumoto 교수님이 작성하신 이 repository 안에는 가스 하이드레이트의 분자 구조를 손쉽게 제작할 수 있는 다양한 모듈들이 있습니다. 다양한 상황에 대비하여 하이드레이트의 구조를 변형하고, 포집되는 기체의 종류도 커스터마이징할 수 있습니다. 때문에 많은 선행 논문들에서 이 repository를 사용하였고, 저 또한 하이드레이트 시뮬레이션을 진행할 때 해당 repository를 유용하게 잘 활용하고 있습니다.

## Installation
GenIce is registered to [PyPI (Python Package Index)](https://pypi.python.org/pypi/GenIce).
Install with pip3.

    pip3 install genice2

## Uninstallation

    pip3 uninstall genice2

$ sudo apt-get install python3-pip -y

$ pip3 install genice2

## Usage

usage: genice2 [-h] [--version] [--rep REP REP REP]
               [--shift SHIFT SHIFT SHIFT] [--dens DENS] [--add_noise percent]
               [--seed SEED] [--format name] [--water model] [--guest 14=me]
               [--Guest 13=me] [--Group 13=bu-:0] [--anion 3=Cl]
               [--cation 3=Na] [--depol DEPOL] [--asis] [--debug] [--quiet]
               [--assess_cages]
               Type

## Clathrate hydrates

For clathrate hydrates, you can prepare the lattice with cages partially occupied by various guest molecules.

* To make a CS1 clathrate hydrate structure of TIP4P water containing CO2 in GROMACS .gro format: (60% of small cages are filled with co2 and 40% are methane)

        genice2 -g 12=co2*0.6+me*0.4 -g 14=co2 --water tip4p CS1 > cs1.gro

* To make a CS2 clathrate hydrate structure of TIP5P water containing THF molecules in the large cage, while only one cage is filled with methane molecule, first, just run `genice` without guest specifications:

        genice2 CS2 > CS2.gro

    The list of cages will be output as follows:

        INFO   Cage types: ['12', '16']
			  INFO   Cage type 12: {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183}
			  INFO   Cage type 16: {136, 137, 138, 139, 140, 141, 142, 143, 16, 17, 18, 19, 20, 21, 22, 23, 160, 161, 162, 163, 164, 165, 166, 167, 40, 41, 42, 43, 44, 45, 46, 47, 184, 185, 186, 187, 188, 189, 190, 191, 64, 65, 66, 67, 68, 69, 70, 71, 88, 89, 90, 91, 92, 93, 94, 95, 112, 113, 114, 115, 116, 117, 118, 119}

(2) 사용 예시

예를 들어, 다음과 같은 형태의 clathrate structure를 만들고 싶다고 가정해봅시다.

• Water type: tip4p model
• Guest type: empty (포집 기체를 안 넣고 싶음)
• Cell replication (unit cell 몇 배의 structure를 만들 것인가?): 2 x 2 x 3 (x축으로 2배, y축으로 2배, z축으로 3배)
• Clathrate type: sI structure
• Output file name: 223sI.gro

코랩 링크에 접속한 후,

그림과 같이 주황색 영역의 “실행” 버튼을 클릭해주세요. 필요한 실행 파일들과 라이브러리들이 자동으로 설치됩니다.

이제 이 셀에서 위에서 언급한 하이드레이트 구조를 만들 수 있도록 합니다. 셋팅을 확인 후, “실행” 버튼을 눌러주세요.

다음 셀에서는 생성한 파일을 미리 Py3Dmol이라는 패키지를 통해 미리 확인해볼 수 있습니다.

마지막 셀을 실행하면, 생성한 구조 파일을 다운로드 해보실 수 있습니다.

다운로드한 223sI.gro 는 visualization software 들을 통해 그 구조를 쉽게 확인해볼 수 있습니다.

대표적인 visualization software들은 다음과 같은 것들이 있습니다.

• OVITO
• VESTA
• PyMOL
• VMD
• Schrodinger

물론, 편한 것을 사용하면 됩니다. 필자는 OVITO 에 익숙하기 때문에 그에 맞춰 설명할 예정입니다.

OVITO를 실행하여 생성된 223sI.gro를 드래그하여 넣어봅니다.

구조가 정상적으로 불러와진 것을 확인할 수 있습니다.

2. GROMACS 구동에 필요한 파일 종류들 알아보기

본격적으로 GROMACS 구동을 진행하기 위해, 필요한 파일 종류들을 모두 알아보겠습니다.

• .itp

  ⇒ Atom, bond, angle, dihedral 등의 정보를 포함하는 파일입니다.

• .atp

  ⇒ Atom types 정보가 기입되어 있는 파일입니다.

• .top

  ⇒ 시스템 내 residue 개수 및 힘장 폴더를 표기하는 파일입니다.

• .log

  ⇒ MD 계산이 완료된 후 간략한 출력 정보를 포함하는 파일입니다.

• .trr

  ⇒ MD 계산 결과가 포함되어 있는 trajectory 파일입니다. MD 분석 시에 대체로 이 파일을 다루게 됩니다.

• .tpr

  ⇒ binary run input file 입니다. MD 계산에 필요한 모든 input과 topology 등이 포함되어 있습니다.

• .mdp

  ⇒ MD 계산을 어떤 식으로 진행할 것인지 사용자가 input 값들을 기입하는 파일입니다.

• .edr

  ⇒ MD 계산 시 출력되는 에너지 정보들이 포함된 파일입니다.

• .cpt

  ⇒ MD 계산을 도중부터 시작할 수 있게 하는 체크포인트 파일입니다.

• .gro

  ⇒ 분자 개수와 각 분자들의 위치 정보들을 포함한 파일입니다.

해당 파일들은 시뮬레이션을 진행하면서 자주 다루게 될 파일들입니다.

당장은 익숙치 않지만, 여러 번 시뮬레이션을 진행하면서 자연스럽게 해당 파일들의 쓰임새를 익히게 되니, 학습과 관련하여 너무 큰 부담을 가질 필요는 없습니다.

3. 간단한 시뮬레이션 진행해보기

(1) .gro 파일 열어보기

1. GenIce 소개 에서 제작한 223sI.gro를 통해 간단한 시뮬레이션을 진행해보겠습니다.

우선 223sI.gro 파일 내부를 간단히 둘러보겠습니다. 해당 파일을 이전 튜토리얼에서 다운로드하였던 텍스트 편집기를 통해 열어봅시다.

Generated by GenIce https://github.com/vitroid/GenIce 
2208
    1ICE     OW    1   1.065   0.370  -0.004
    1ICE    HW1    2   1.037   0.319   0.071
    1ICE    HW2    3   1.010   0.448  -0.003
    1ICE     MW    4   1.055   0.374   0.006
    2ICE     OW    5   0.370   0.004   0.135
    2ICE    HW1    6   0.448   0.003   0.191
    2ICE    HW2    7   0.319  -0.071   0.163
    2ICE     MW    8   0.374  -0.006   0.146
    3ICE     OW    9   0.005   0.300   0.605
.
.
.
	552ICE     OW 2205   2.030   1.205   3.466
  552ICE    HW1 2206   2.082   1.129   3.438
  552ICE    HW2 2207   1.953   1.203   3.410
  552ICE     MW 2208   2.027   1.195   3.455
    2.4006977272287195 2.4006977272287195 3.6010465908430795
#

• 2208은 해당 시스템 안에 포함된 원자의 총 갯수입니다.
• ICE는 GenIce 모듈이 임의로 작명한 물 분자의 residue 이름 입니다.
• ICE앞의 숫자는 residue 번호입니다.
• OW, HW1, HW2, MW은 각각 TIP4P 물 분자들을 구성하고 있는 원자들의 종류입니다.
• 오른쪽의 3열 숫자들은 각각 해당 원자들의 x, y, z 좌표를 나타냅니다.
• 2.4006977272287195 2.4006977272287195 3.6010465908430795은 이 시스템의 boundary를 나타냅니다. 즉, 2.4 x 2.4 x 3.6 $n{m}^3$의 크기를 가진 상자라는 의미입니다.

(2) 힘장 (Force fields) 이란?

이제 힘장에 대해서 이야기 해보도록 하겠습니다.

위키피디아에 따르면, 힘장은

A force field is a computational method that is used to estimate the forces between atoms within molecules and also between molecules. More precisely, the force field refers to the functional form and parameter sets used to calculate the potential energy of a system of atoms or coarse-grained particles in molecular mechanics, molecular dynamics, or Monte Carlo simulations.

라고 언급되어 있습니다. 즉, 힘장은 MD 계산을 위해 설정해야 할 파라미터 값들입니다.

그렇다면, 어떤 파라미터 값들을 설정해야 MD 시뮬레이션을 돌릴 수 있을까요?

Bonds, Angles, Non-bonded angles, 그리고 charges가 해당 파라미터 값들입니다.

(3) 힘장 폴더 로컬 디렉토리로 불러오기

다행히도, 이 모든 파라미터 값을 손으로 작성할 필요는 없습니다.

GROMACS가 설치된 폴더에는, 다양한 종류의 ff (force fields) 폴더들이 있습니다. 이를 템플릿 삼아 조금씩 변형하면 됩니다.

우선 이 문서에서는 OPLS 힘장을 활용해보겠습니다.

다음 커맨드를 입력하여, GROMACS 폴더에서 로컬 디렉토리로 OPLS 힘장 폴더를 불러옵시다.

$ cp -r /usr/local/gromacs/share/gromacs/top/oplsaa.ff ./oplsaa.ff

현재 작업 중인 로컬 디렉토리에 oplsaa.ff/라는 이름의 폴더가 생성되는 것을 확인할 수 있습니다.

이 폴더 안에는 다음과 같은 파일들이 있음을 확인할 수 있습니다.

spc, tip3p, tip4p 등 다양한 종류의 water model들이 있습니다.

(4) 힘장 파일 (.itp) 및 분자 파일 (.gro) 커스터마이징

불행히도, 위에서 언급한 .gro, oplsaa.ff/ 만으로는 계산이 불가능합니다.

그 이유는, Genice 소스 코드에서 생성해 낸 물 분자 정보와 oplsaa.ff/ 내의 tip4p.itp 의 원자 이름 및 배열 순서가 완전히 일치해야 GROMACS에서 이를 인식하기 때문입니다.

이를 해결하기 위해서는, 다음 절차를 진행하세요.

1. 코랩 링크에 다시 접속하여, 이번에는 아래 그림과 같이 water_model을 tip4p가 아닌, cttip4p로 바꿔서 하이드레이트 구조를 만들어주세요.

   

2. 로컬 컴퓨터로 돌아와, 앞서 복사해놓았던 oplsaa.ff/에 진입하여 tip4pew.itp를 다음과 같이 수정합니다. .itp 파일을 텍스트 편집기 등을 통해 열어 내용을 수정해야 합니다.

   [ atoms ]
   ; id  at type     res nr  res name  at name  cg nr  charge    mass
     1   OW_tip4pew  1       SOL       OW       1       0        16.00000
     2   HW_tip4pew  1       SOL       HW1      1       0.52422   1.00800
     3   HW_tip4pew  1       SOL       HW2      1       0.52422   1.00800
     4   MW          1       SOL       MW       1      -1.04844   0.00000

   [ atoms ]
   ; id  at type     res nr  res name  at name  cg nr  charge    mass
     1   OW_tip4pew  1       SOL       O       1       0        16.00000
     2   HW_tip4pew  1       SOL       H      1       0.52422   1.00800
     3   HW_tip4pew  1       SOL       H      1       0.52422   1.00800
     4   MW          1       SOL       M       1      -1.04844   0.00000

이렇게 하면, oplsaa.ff/ 폴더 내에 있는 tip4p.itp와 223sI.gro 파일이 GROMACS에서 호환됩니다.

(5) .grompp 키워드를 이용하여 GROMACS 구동 준비하기

이제 본격적으로 GROMACS를 구동해보겠습니다.

gmx grompp를 이용하면 .mdp파일을 통해 입력한 다양한 셋팅 값들을 준비한 .gro파일과 병합시켜 .tpr 파일을 제작해 GROMACS를 구동할 준비를 완료할 수 있습니다.

(보다 상세한 사용법은 gmx grompp -h 를 치면 터미널 내에서 사용 방법을 볼 수 있습니다)

이 튜토리얼에서는 미리 작성이 완료된 .mdp 파일들을 이용하여 gmx grompp 명령어를 사용합니다.

먼저, 아래 파일들을 받아 로컬 디렉토리에 다운로드 해주세요. 이 튜토리얼에서 진행해 볼 시뮬레이션 셋팅이 작성된 파일들입니다.

em.mdp 1.0KB

npt.mdp 2.2KB

nvt.mdp 2.1KB

보통 논문에서 computational methods 라고 적힌 섹션에 언급되어 있는 시뮬레이션 방법들은 대부분 이 .mdp 파일들을 이용하여 셋팅하게 됩니다. 대부분의 셋팅은 대부분 건들지 않고 계산하는 것이 일반적입니다. 저희가 주로 건드리게 될 항목들은 다음과 같습니다.

• rcoulomb - Short-range electrostatic cut-off
• rvdw - Short-range Van der Waals cut-off
• nsteps - 시뮬레이션 길이를 결정
• dt - 몇 timestep으로 계산할 것인지?
• Output control
nstxout save coordinates every xxx steps
  nstvout save velocities every xxx steps
  nstenergy save energies every xxx steps
  nstlog update log file every xxx steps
• tc-grps - 어떤 residue의 온도를 조절할 것인가?
• ref_t - 온도 설정은? (단위: Kelvin)
• ref_p - 압력 설정은? (단위: bar)

그 외에도 많은 옵션들이 있으므로, 상세한 내용은 GROMACS manual을 참고하시면 됩니다.

이제 마지막으로, 아래 파일도 다운로드 받아주세요.

topol.top 0.2KB

.top 파일은 우리가 시뮬레이션 할 시스템의 stoichiometry 정보와 힘장 정보를 포함한 파일입니다.

.gro 파일 안에는 552개의 SOL 분자들이 있으므로 다음과 같이 입력해주면 됩니다.

[ molecules ]
; Compound	#mols
SOL			552

이제 gmx grompp 를 실행합니다.

$ gmx grompp -f em.mdp -c ct223sI.gro -p topol.top -o em.tpr 

아래와 같이 메시지가 출력됩니다.

Setting the LD random seed to -2013923346

Generated 330891 of the 330891 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5

Generated 330891 of the 330891 1-4 parameter combinations

Excluding 2 bonded neighbours molecule type 'SOL'

Cleaning up constraints and constant bonded interactions with virtual sites
Analysing residue names:
There are:   552      Water residues
Number of degrees of freedom in T-Coupling group rest is 3309.00
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 24x24x32, spacing 0.100 0.100 0.113

Estimate for the relative computational load of the PME mesh part: 0.24

This run will generate roughly 0 Mb of data

디렉토리를 확인해보면, em.tpr이라는 파일이 새롭게 생겼습니다. 저희가 만든 시스템의 energy minimization을 구동할 준비가 완료되었다는 뜻입니다.

이 .tpr 파일이 있으면 이제 시뮬레이션을 시작할 수 있습니다. 다음 명령어를 입력합니다.

$ gmx mdrun -v -deffnm em

Steepest Descents:
   Tolerance (Fmax)   =  1.00000e+03
   Number of steps    =        50000
Step=    0, Dmax= 1.0e-03 nm, Epot= -3.21428e+04 Fmax= 3.83272e+02, atom= 1621

writing lowest energy coordinates.

Steepest Descents converged to Fmax < 1000 in 1 steps
Potential Energy  = -3.2142805e+04
Maximum force     =  3.8327237e+02 on atom 1621
Norm of force     =  1.2046063e+02

GROMACS reminds you: "You see it through a charmed medium: you can not discern that the gilding is slime and the silk draperies cobwebs; that the marble is sordid slate, and the polished woods mere refuse chips and scale bark." (Mr. Rochester in Jane Eyre by Charlotte Bronte)

한 스텝만에 energy minimization이 완료되었습니다.

이제 디렉토리를 확인해보면, em.gro, em.trr, em.edr, em.log 파일들이 생성되었습니다.

이 파일들은 모두 시뮬레이션 이후에 출력되는 output 파일들입니다.

곧이어, nvt.mdp를 이용해 NVT relaxation을 진행해보도록 하겠습니다. 아래 커맨드를 입력해주세요.

$ gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr

다음과 같은 메시지가 출력됩니다.

.
.
.
Estimate for the relative computational load of the PME mesh part: 0.24

This run will generate roughly 3 Mb of data

디렉토리에는 nvt.tpr가 생성되어 있습니다.

아래 커맨드를 입력해 NVT relaxation을 시작하겠습니다.

$ gmx mdrun -v -deffnm nvt

아래와 같이 계산이 진행되고,

starting mdrun 'Tutorial1'
25000 steps,     50.0 ps.
step 17100, remaining wall clock time:     6 s

Writing final coordinates.
step 25000, remaining wall clock time:     0 s
               Core t (s)   Wall t (s)        (%)
       Time:      320.688       26.725     1200.0
                 (ns/day)    (hour/ns)
Performance:      161.655        0.148

계산이 종료됩니다.

디렉토리에는 nvt.log, nvt.gro, nvt.trr, nvt.log가 새롭게 생성되었습니다.

이제 relaxation된 시스템과 npt.mdp를 이용해 1 ns 길이의 NPT 시뮬레이션을 진행해보겠습니다.

다음 명령어를 입력해주세요.

$ gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr 

디렉토리에 npt.tpr 파일이 형성됩니다.

다음 명령어를 입력하여 production run (본격적으로 진행할 시뮬레이션을 지칭)을 진행하겠습니다.

$ gmx mdrun -v -deffnm npt

10분 정도 시간이 지나니, 시뮬레이션이 종료되었습니다.

(6) 계산 결과 확인해보기

a. Trajectory 시각화

우선, 시뮬레이션이 제대로 이루어졌는지 확인하기 위해 visualization (시각화)를 진행해야 합니다.

다음처럼 bash script 파일 하나를 만듭시다.

$ vi trjconv.sh

esc - :wq를 입력하고, 텍스트 편집기로 trjconv.sh를 수정합니다.

#!/bin/bash

gmx trjconv -f npt.trr -s npt.tpr -o npt_trjconv_water.gro << EOF
1
EOF

exit

gmx trjconv는 .trr 파일을 visualization이 가능한 .gro 파일로 그 형태를 바꿔주는 GROMACS 내장 모듈입니다.

EOF는 End of File의 약자로, EOF가 다시 작성될 때까지 커맨드를 지속적으로 입력한다는 의미입니다.

여기서는 gmx trjconv -f npt.trr -s npt.tpr -o npt_trjconv_water.gro 명령어를 입력하고, 곧이어 1 번 옵션을 선택하라는 의미입니다.

그 다음, 터미널에서 다음과 같이 bash script를 실행시킵니다.

$ ./trjconv.sh

약 15MB의 npt_trjconv_water.gro 파일이 형성됩니다.

OVITO, VMD와 같은 visualization software에 해당 파일을 넣어보세요.

250 K, 100 bar에서는 속이 빈 sI hydrate의 구조가 녹지 않고 유지되는 것을 육안으로 확인할 수 있습니다.

b. Potential energy 확인하기

이제 수치 자료를 통해 하이드레이트가 녹지 않았다는 것을 증명해보겠습니다.

$ vi poteng.sh

esc - wq:

아래와 같이 내용을 수정합니다.

#!/bin/bash

gmx energy -f npt.edr -o poteng.xvg << EOF
5
EOF

exit

⇒ 여기서는 5번이 potential energy이지만, 다른 시뮬레이션의 경우에는 5번이 아닐 수도 있으므로 주의하세요.

Potential energy를 npt.edr로부터 도출하는 bash script를 완성하였습니다. 이제 이를 실행합니다.

$ ./poteng.sh

로컬 디렉토리에 poteng.xvg 파일이 형성됩니다.

이를 Excel, sigmaplot 등으로도 확인해볼 수 있지만, 저희는 리눅스에서 이를 진행해보겠습니다.

xmgrace라는 툴을 이용해서 .xvg 파일의 결과를 간단히 확인해볼 수 있습니다.

xmgrace는 다음과 같이 설치할 수 있습니다. (왼쪽 화살표를 눌러 상세한 방법을 확인하세요)

sudo systemd-machine-id-setup
sudo dbus-uuidgen --ensure

cat /etc/machine-id 

sudo apt-get install x11-apps xfonts-base xfonts-100dpi xfonts-75dpi xfonts-cyrillic

export DISPLAY=:0
source ~/.bashrc

xeyes 

$ pip3 install grace

$ xmgrace 

설치가 완료되면, 터미널에 다음과 같이 명령어를 입력하세요.

$ xmgrace poteng.xvg

XMGRACE를 통해 1 ns (= 1000 ps)의 시간 동안 potential energy가 -29300 ~ -29800 kJ/mol 사이의 값에서 잘 유지되었음을 확인할 수 있습니다.

c. Radial distribution function (RDF) 그려보기

H. Jiang et al., Molecular Dynamics Simulations of Methane Hydrate Using Polarizable Force Fields 논문에 있는 radial distribution function은 다음과 같습니다.

즉, 하이드레이트의 구조가 유지된다면 물 분자의 산소와 다른 물 분자의 산소 사이의 RDF는 위 그림과 유사해야합니다.

과연 그럴지 실제로 분석을 진행해보도록 하겠습니다.

우선, 물 분자의 산소 원자만을 인덱싱해야 합니다. 인덱싱을 할 때에는 gmx make_ndx 모듈을 사용합니다.

먼저, 다음처럼 Ow_indexing.sh 파일을 만들어줍시다.

$ vi Ow_indexing.sh

esc - wq:

입력 후 내부를 다음 내용으로 수정합니다.

#!/bin/bash

gmx make_ndx -f npt.gro -o Ow_index.ndx << EOF
a O 
q
EOF

그 다음 해당 스크립팅을 실행합니다.

$ ./Ow_indexing.sh 

로컬 디렉토리에 Ow_index.ndx 파일이 생겼습니다.

이제 gmx rdf 모듈을 활용하여 0 ~ 200 ps, 그리고 800 ~ 1000 ps 에서의 $O_w-O_w$ RDF ($g_{O_w-O_w}(r)$) 를 비교해보겠습니다.

(프로파일이 비슷하다면 하이드레이트 구조가 잘 유지된 것, 그렇지 않다면 잘 유지되지 않은 것)

rdf_making.sh 파일을 만듭니다.

$ vi rdf_making.sh

다음과 같이 rdf_making.sh의 내용을 수정해줍니다.

mkdir rdf_comparison/ &&

gmx rdf -f npt.trr -s npt.tpr -o ow-ow_0-200ps.xvg -n Ow_index.ndx -tu ps -b 0 -e 200 << EOF
3
3
EOF

gmx rdf -f npt.trr -s npt.tpr -o ow-ow_800-1000ps.xvg -n Ow_index.ndx -tu ps -b 800 -e 1000 << EOF
3
3
EOF

mv ow-ow_0-200ps.xvg rdf_comparison/ &&
mv ow-ow_800-1000ps.xvg rdf_comparison/ 

3 번 옵션은 Ow_index.ndx 파일로 인덱싱 한 산소 원자들 (O)을 의미합니다.

즉, O와 O 사이의 RDF를 npt.trr, npt.tpr 을 통해 도출한다는 명령입니다.

-tu 옵션은 time unit을 의미합니다. 이 시뮬레이션에서는 ps를 단위로 놓았습니다.

-b 옵션은 begin time, 즉 RDF 분석을 시작할 시간 지점,-e 옵션은 end time, RDF 분석을 끝낼 시간 지점을 각각 의미합니다.

이제 스크립트를 실행하여 rdf 분석을 진행합니다.

$ ./rdf_making.sh

로컬 디렉토리에 rdf_comparison/라는 이름의 폴더가 생성되었고, 이 안에 ow-ow_0-200ps.xvg, ow-ow_800-1000ps.xvg 이렇게 두 파일이 만들어져 있습니다. 이제 다음 커맨드를 이용하여 두 RDF 그래프를 비교해봅니다.

$ cd rdf_comparison

$ xmgrace *.xvg

*(asterisk) 는 리눅스에서 임의의 이름을 지칭할 때 유용합니다. 여기서는 rdf_comparison/ 폴더 내의 두 파일을 지칭합니다.

xmgrace를 활용하면 다음과 같은 그림이 그려집니다.

이제 레퍼런스와의 비교를 위해, 그래프의 크기와 스케일 조정을 진행해주겠습니다. xmgrace에서의 그래프 수정은 간단합니다.

Axis properties... 를 누르고,

이와 같은 창에서 그래프 속성을 이리저리 변경해주면 됩니다.

• 레퍼런스와 사용 힘장 및 unit cell 크기 등이 다르기 때문에 프로파일이 조금 다를 수는 있지만, 큰 경향은 대체로 일치하는 것을 알 수 있습니다.
• 검은 색 (0 ~ 200 ps) 프로필과 붉은 색 (800 ~ 1000 ps)의 프로파일이 거의 일치하므로, 하이드레이트 구조가 변화하지 않고 (여기서는 녹지 않고) 잘 유지되었음을 확인할 수 있습니다.

◻️