NAMD, VMD 튜토리얼

1 개요[ | ]

Tutorials on NAMD with VMD
NAMD, VMD 튜토리얼

by One-Sun Lee, currently at QEERI

2 시작하기[ | ]

필요한 파일들:

NAMD	http://www.ks.uiuc.edu/Research/namd/
VMD	http://www.ks.uiuc.edu/Research/vmd/

NAMD는 실제로 계산을 하는 엔진이고, VMD는 이것에 필요한 인풋 파일을 만들고 계산과정을 중간중간에 점검하며 콘트롤하고 결과를 분석하는 프로그램이다. 윈도우용과 리눅스를 비롯한 각종 유닉스용 모두 다운 받을 수 있다.

force field는 CHARMm을 사용하겠다. 사실은 CHARMm, AMBER, GROMOS가 가능하지만 일단은 CHARMm으로 시작하자. CHARMm force field는 http://www.pharmacy.umaryland.edu/~alex/research.html 에 가면 구할 수 있다.

첫번째 예제는 BPTI분자를 pdb에서 가져와서 수소를 붙인 후, N-terminal에는 ACE group을 붙이고 C-terminal에는 CT2를 붙여 보는 것이다. 이 예제는 NAMD user guide에 있는 것을 약간 고친 것이다. directory 구조는 편의상 current directory가 bpti이고 그 아래에 toppar과 output directory가 있는것으로 하자.

1. toppar에 force field file들을 넣어두자. 여기에 필요한 파일은 top_all22_prot.inp 과par_all22_prot.inp 이다. 새로운 버젼을 사용해도 된다. (사실 지금은 27도 나와 있다.)

2. 그런후 http://www.rcsb.org/pdb/ 에 가서 분자를 가져오자.

찾기에서 6pti라고 친후 Bovine Pancreatic Trypsin Inhibitor (BPTI, Crystal Form III) 라는 것을 가져오면 된다. (맨 아래 나온다.) 여기에는 수소의 위치는 나와 있지 않다. 이걸 가져와서 파일 내용을 들여다 보면, 이 안에 단백질의 구조만 있는 것이 아니라 약간의 물 분자도 섞여 있는데 NAMD 파일을 만들려면 일단 이 두개를 분리해야 한다. 분리하는 방법은

grep -v '^HETATM' 6PTI.pdb > output/6PTI_protein.pdb
grep 'HOH' 6PTI.pdb > output/6PTI_water.pdb

3. NAMD에서 계산을 하려면 .psf와 .pdb 파일이 필요한데 .psf psfgenerator나 VMD를 사용해서 만들수 있다. 여기서는 VMD를 사용해서 만드는 것을 보자. (사실은 이 두개의 방법은 동일하다. 왜냐면 VMD안에 psfgenerator가 내장되어 있어서 사실은 psfgenerator를 사용하기 때문이다.)

VMD를 설치한 후 실행하면 3개의 창이 뜨는데 그 중 command창이 있다. directory를 이미 만들어 놓은 bpti로 옮긴다. 그런 후

topology toppar/top_all22_prot.inp

라고 치면 topology file을 읽어온다. 그 후에

segment BPTI {
pdb output/6PTI_protein.pdb
first ACE
last CT2
}

라고 치면 BPTI라는 segment를 만들고 이 안에 단백질의 sequence정보를 읽어 들인 후 N-terminal과 C-Terminal에 각각 ACE와 CT2를 patch한다. 이 ACE 와 CT2는 Patch RESidue (PRES)의 일종인데 topology file을 열어보면 그 안에 정의되어 있는 것을 볼 수 있다. PRES중 다른 것을 골라도 된다.

4. sulfide bond를 아래와 같이 patch한다.

patch DISU BPTI:5 BPTI:55
patch DISU BPTI:14 BPTI:38
patch DISU BPTI:30 BPTI:51

그리고 ILE residue에서는 CD의 이름이 CD1이므로 alias가 필요하다.

alias atom ILE CD1 CD

5. 이제 원자의 좌표를 읽어보자.

coordpdb output/6PTI_protein.pdb BPTI

라고 치면 아까 만들어 두었던 BPTI라는 segment에 6PTI_protein.pdb의 좌표를 읽어온다.

6. 이제 새로운 psf 파일과 pdb파일을 만들자.

writepsf output/bpti.psf

라고 치면 output및에 bpti.psf라는 psf 파일이 만들어지고, 다시

guesscoord

이라고 치면 새로 만들어 준 구조에서 수소원자들의 위치를 적당히 계산한다. (왜냐하면 아까 coordpdb에서 읽어온 좌표에는 수소의 위치가 없으므로.) 그리고

writepdb output/bpti.pdb

라고 치면 수소원자를 붙여서 valence를 맞춘 pdb 파일이 생성된다.

3 단백질의 간단한 dynamics 계산 및 분석 (1)[ | ]

만들어 놓은 bpti.psf와 bpti.pdb를 가지고 CHARMm force field를 이용해 간단한 dynamics simulation을 해보자.

여기에 또하나 필요한 것은 configuration file인데, 이것이 실제 명령어 모음 파일이다. 즉, 인풋 파일은 이름이 뭐고, 아웃풋은 이름이 뭐며, 어디서 쓸것인가, 얼마만큼 minimize하고 dynamics를 할 것인가, 원자의 전하는 어떻게 계산할 것이가 등등 모든 계산에 관련된 명령어가 여기에 모두 모아져 있다.

configure file의 간단한 예를 한 번 보자. (이 파일의 이름은 아무렇게나 써도 된다. 즉, 확장자 걱정은 안해도 된다.)

# NAMD configuration file for BPTI
# molecular system
structure output/bpti.psf
# force field
paratypecharmm on
parameters toppar/par_all22_prot.inp
exclude scaled1-4
1-4scaling 1.0
# approximations
switching on
switchdist 8
cutoff 12
pairlistdist 13.5
margin 0
stepspercycle 20
#integrator
timestep 1.0
#output
outputenergies 10
outputtiming 100
binaryoutput no

# molecular system
coordinates output/bpti.pdb
#output
outputname output/bpti_out1
dcdfreq 1000
#protocol
temperature 0
reassignFreq 1000
reassignTemp 25
reassignIncr 25
reassignHold 300
#script
minimize 1000
run 20000

각각의 부분에 대한 자세한 설명은 나중에 하고 일단 계산을 시작해 보자. 계산은 윈도우에서 도스창을 띄우고 해도 되고, 리눅스나 기타 유닉스에서 해도 된다. 사실 NAMD는 parallel 계산이 강력하지만 지금은 single 계산으로만 해보자. 명령은

$path/charmrun ++local +pN $path/namd2 configure_file_name > output_file_name &

으로 하면된다. 여기서 N은 사용하고자하는 processor의 갯수이다. 따라서 만일 2개의 processor를 사용하고 싶으면 +p2라고 써 주라는 뜻이다. 그러나 어떤 오퍼레이팅 시스템인가에 따라서 또는 qsub등의 관리자를 사용하는것에 따라서 약간씩 옵션이 다르게 사용되므로 자세한것은 메뉴얼에서 찿아보자. 다만 이때 path에 주의하자. 만일 계산하는 디렉토리가 configure file이 있는 곳이면 (이렇게 하는 것이 좋다. 왜냐하면 charmrun이 있는 곳에서 시작하면 여기다 결과를 쓰게 된다.) 미리 .cshrc 파일등에 charmrun과 namd2의 path를 지정해야 한다.

계산이 끝나면 .dcd file이 생성되는데 이것을 보려면 VMD를 사용해야 한다. VMD를 실행하고 옥색 (파란색?)창에서 load를 눌러보면 primary file type과 secondary file type이 있는데 앞에는 h2c.psf를, 뒤에는 h2c_out1.dcd를 선택한다. 그러면 일종의 movie를 볼 수 있는데 여기에 속도 조절하는 것도 있으니 적절히 선택해 보면 된다.

4 단백질의 간단한 dynamics 계산 및 분석 (2)[ | ]

단백질의 Molecular Dynamics 계산 후 가장 많이 쓰이는 분석 방법 중 하나가 Root Mean Square Deviation (RMSD) 를 보는 것이다. 앞에서 설명한 Molecular Dynamics계산후 .DCD file을 얻게 되는데 이 파일을 이용하여 RMSD를 구해보자.

1. 우선 VMD에서 계산이 끝난 결과를 불러온다. 이때 primary구조에선 psf파일을 불러오고, secondary구조에선 dcd 파일을 불러온다.

2. 그 후 다음과 같은 스크립트를 VMD command 화면에서 실행한다.

set reference [atomselect top "protein" frame 0]
set compare [atomselect top "protein"]
set num_steps [molinfo top get numframes]
set out [open test_md.dat w]
for {set frame 0} {$frame < $num_steps} {incr frame} {
	$compare frame $frame
	$compare move [measure fit $compare $reference]
	puts $out "$frame  [measure rmsd $compare $reference]"
}
close $out

이것은 맨 처음 frame을 기본으로 해서 그 이후의 프레임들의 RMSD를 구하는것이다. 이렇게 하면 test_md.dat라는 파일을 만들고 거기에 각각의 프레임별로 RMSD를 계산해준다. 결과는

0     0
1     1.328801274
2     1.610798717
3     2.366841793
4     2.613398314
5     2.49978447
6     2.593734741
7     2.644852638
8     2.791284084
9     2.591056108
10     2.742631674
... (생략)

와 같은 모양의 파일을 얻게 된다. 첫번째 컬럼은 프레임넘버이고 두번째 컬럼이 각각의 RMSD (단위는 거리) 값이다.

Upload:rmsd.gif

5 두 개 이상의 분자를 계산할 때[ | ]

두 개 이상의 분자를 계산하는 것은 한개의 분자를 계산할 때와 기본적인 명령어는 동일하지만, namd 에서 pdb와 psf 파일을 가지고 계산을 할 때에는 고려해 줘야 할 문제점이 있다.

보통 단백질을 계산할 때 엑스레이 결정으로부터 시작하는데, 이 엑스레이 결정 pdb 파일에는 수소원자에 대한 정보가 없다. 따라서 이 엑스레이 결정 구조 정보에서 psf 파일과 수소원자를 붙여서 옥텟룰을 맞춘 pdb 파일을 다시 만들어야 한다.

그런데 이런 경우를 생각해보자: 크리스탈 구조가 두개의 분자로 이루어져 있는데 (편의상 분자 A와 B라고 하자), 분자 A는 residue number가 1부터 160까지있고, B는 500부터 660까지 있고 가운데는 없다 (이런 경우는 아주 흔하다). 또 우리는 1번과 500번에 ACE를 붙이고 160번과 660번에 CT3를 붙이고 싶어한다고 가정하자.

그냥 앞에서 설명한대로 VMD에서 psf 파일을 만들면 VMD는 이 두 분자를 하나로 인식해서 residue 1번에 ACE를 붙이고 660번에 CT3를 붙인후 residue 160과 500을 연결시켜 psf에 쓴다.

이 경우에는 일단 pdb파일을 에디터를 이용해 두 개의 파일로 자른다. (지금의 예에서는 1 ~ 160까지 하나의 파일, 그리고 500 ~ 660까지 또 하나) 그리고 나서 한 개만 먼저 VMD에서 불러서 ACE와 CT3를 patch한다. VMD를 끄지 말고 계속해서 (끄지 말아야한다 !) 다음 번 파일을 불러와서 다시 ACE와 CT3를 patch한다. 그런 다음 writepsf와 writepdb을 사용하여 psf와 pdb 파일을 만든다. 그러면 각각 분자에 ACE와 CT3 terminal을 붙인 결과를 얻을 수 있다.

사실 이 방법은 좀 이상하기는 하다. 뭔가 더 쉬운 방법이 있을 듯도 싶은데. 이 방법은 메뉴얼 어디에도 없다. 오늘 오후 내내 혼자 이리저리 해 본 결과다. 누가 더 쉬운 방법 있음 여기에 알려 줘라. (아마도 pdb 파일에 그냥 정보를 주면 될 것 같은데 어떻게 써야하는지 모르겠다. 일단은 그냥 쓰자.)

6 Charge fitting[ | ]

원자의 전하(atomic charge)는 분자간 힘에서 중요한 부분을 차지한다. 대부분의 분자의 모양은 원자의 전하에 가장 크게 영향을 받으므로 force field를 사용하는 계산에서 원자의 전하값은 되도록이면 정확하게 지정하여야 한다. 그래서 NAMD/VMD와는 직접적으로 관계가 없지만 원자의 전하 계산에 관한 내용을 언급하고자 한다.

NAMD/VMD에서는 원자의 전자를 직접 계산할 수 없기 때문에, 다른 프로그램을 사용하여 원자의 전하를 계산한 후 직접 집어 넣어 줘야 하는 경우가 종종 발생한다. 단백질의 경우는 미리 atomic charge가 지정되어 있지만 그 이외의 경우는 값이 없기 때문이다.

atomic charge는 다른 물리량과는 달리 직접 실험적으로 정해줄 방법이 없으므로 여러가지 이론적인 방법이 개발되어 있다. 그 중 가장 많이 사용되는 것이 Gaussian이나 GAMESS를 사용하여 Eletrostatic Potential 을 계산한 후 atomic charge를 이것에 맞도록 fitting하는 방법이다. 여기서는 Kollman group에서 개발한 RESP (Restricted ElectroStatic Potential)을 설명한다.

RESP는 Kollman group에서 개발한 AMBER안에 들어있는 모듈 중 하나인데, AMBER는 상용 프로그램이지만 RESP는 따로 떼어서 공짜로 나누어 준다. Kollman group의 홈페이지는 http://www.amber.ucsf.edu/amber/amber.html 이다. 여기서는 가장 간단한 예로 pyridine의 atomic charge를 계산해보자. pyridine의 모양과 각각의 원자의 번호는 그림에 있다. 여기서 번호는 원자번호가 아니라 편의상 붙인 이름이다.

Upload:charge.gif

위에서 보듯이 탄소 2번과 6번은 대칭의 위치에 있으므로 서로 전하가 동일하여야 한다. 마찬가지 이유로 탄소 3과 5, 수소 8과 11, 수소 9와 10도 전하가 동일하여야 한다. 그런데 일반적으로 그냥 Gaussian 등의 프로그램을 사용하여 atomic charge를 계산하면 이들의 전하가 동일하지 않게 나오는 경우가 있다. 이것은 주위의 환경이 약간만 달라도 계산이 영향을 받기 때문이다. 물론 위의 그림에서는 주위에 다른 분자 또는 원자가 없으므로 이런 경우에는 동일한 전하를 얻게 되지만...여기서는 하여간 피리딘이 가장 간단한 분자 가운데 하나이므로 예제로 택했다.

일단 Gaussian geometry optimization을 하자. 인풋파일은 다음과 같다.

# HF/6-31G* Opt
RESP test
0 1
C
C,1,R2
N,1,R3,2,A3
... (생략)

geometry optimization이 끝나면 그 결과 구조를 가지고 다시 다음과 같은 계산을 한다.

# HF/6-31G* pop=MK iop(6/33=2)
RESP pyridine test
0 1
... (생략)

여기서 반드시 pop=MK iop(6/33=2)의 옵션을 집어넣어야 한다. 그리고 geometry optimization에서과 동일한 basis set을 사용하여야 한다. 그런 후 RESP가 있는 디렉토리에 가보면 esp.sh라는 스크립트가 있다. 이 스크립트와 gaussian결과를 이용하여 다음과 같이 실행한다. (RESP는 Linux를 포함하여 UNIX계열에서만 실행된다.)

./esp.sh gaussian_output_file

그러면 esp.dat라는 파일이 생성된다. 이제 RESP의 첫 번째 인풋 파일을 만들어 보자. 첫 번째 인풋파일 (resp1.in) 은 아래와 같다.

input pyridine HF/6-31G*
&cntrl nmol=1, ihfree=1
&end
1.0
tol_resp1
	0   11
	6    0
	6    0
	6    0
	7    0
	6    0
	6    0
	1    0
	1    0
	1    0
	1    0
	1    0

아시다시피 인풋 파일을 작성할 때에는 줄을 잘 맞춰야 한다.

첫번째 줄은 그냥 코멘트이고 2 ~ 4번째 줄은 그냥 그래도 넣어주면 된다. 5번째 줄도 코멘트이다.

그 다음에 나오는 0 11은 전체 전하가 0 이고 원자의 갯수가 11개라는 의미이다. 그 다음 줄부터는 각각의 원자에 대한 정보인데 한 줄이 원자 한 개에 대한 정보이다. 처음에 나온 6, 7, 또는 1이라는 숫자는 atomic number즉 탄소, 질소, 수소를 의미하고 그 다음에 붙어 있는 0은 charge fitting을 어떻게 해 줄 것인가에 대한 정보이다. 첫 번째 input에서는 대체로 0으로 놓고 계산을 실행해도 큰 문제는 없다. 그런 다음 다시 RESP가 있는 디렉토리로 돌아와서 다음과 같이 실행한다.

./resp -i resp1.in -o resp1out -p punch1resp -e esp.dat -t resp1.chg

여기서 resp1.in은 위에서 만든 파일이고 esp.dat또한 아까 위에서 얻은 파일이다. 그 이외의 파일들은 (resp1out, punch1resp, resp1.chg) 위 명령을 실행하면 만들어진다. 그런 다음 두번째 인풋 (resp2.in) 을 만들어보자.

input pyridine HF/6-31G*
&cntrl nmol=1, ihfree=1, iqopt=2, qwt=0.001
&end
1.0
tol_resp1
	0   11
	6  -99
	6    0
	6    0
	7  -99
	6    3
	6    2
	1  -99
	1    0
	1    0
	1    9
	1    8

여기에서는 iqopt=2를 반드시 넣어주어야 한다. 위 인풋에서 일번 탄소에 해당하는 곳에 6 -99라고 써 있는 것은 이 탄소에 대한 전하는 첫번째 계산에서 사용한 것을 바꾸지 않겠다는 뜻이다. 두번째 탄소에 해당하는 곳에 있는 6 0은 이 탄소는 전하계산을 다시 하겠다는 뜻이고, 6번 탄소에 해당하는 곳에 있는 6 2는 6번 탄소의 전하는 2번탄소와 동일하게 만들라는 의미이다. 마찬가지로 11번 수소에 해당하는곳에 있는 1 8는 11번 수소는 8번 수소와 전하게 동일하게 만들라는 의미이다. 나머지도 마찬가지로 생각하면 된다.

이 파일은 만들고 난뒤 다음과 같이 실행한다.

./resp -i resp2.in -o resp2out -e esp.dat -t resp2.chg -p punch2resp -q resp1.chg

여기서 사용된 resp1.chg는 위 첫번째 인풋에서 만들어진 파일이다. 위처럼 실행하면 punch2resp가 생성되는데 이것이 우리가 원하는 파일이다. 이 안에 새롭게 fitting된 전하에 대한 정보가 있다.

7 분자를 z-축에 정렬하고 z-축으로만 움직이는 MD하기[ | ]

lipid 사이에 들어 있는 단백질 channel 등을 만들 때, 종종 이 단백질을 z-축에 정렬해야 하는 일이 생긴다. NAMD의 경우, 이러한 종류의 모의 실험을 위하여 package가 만들어진 것이 있다. 이것은 배포되는 VMD에는 들어있지 않고 VMD website에서 따로 받아서 plugin 해야한다. 필요한 파일은 http://www.ks.uiuc.edu/Research/vmd/script_library/scripts/orient/ 에서 받을수 있다.

우선, 여기에서 orient.tar.gz 와 la101psx.tar.gz 2개의 파일을 받아서 압축을 푼다. 압축을 풀어서 디렉토리를 $VMD_HOME/plugins/noarch/tcl 아래에 만든다. 그리고 VMD를 실행시켜 원하는 분자를 연다. 그런 다음 VMD command 창으로 가서

package require La
package require Orient

라고 입력하면 각각의 모듈의 버젼을 보여주면서 plugin 되는 것을 볼 수 있다. 그 다음 아래 명령을 실행한다.

namespace import Orient::orient
set sel [atomselect top "all"]
set I [draw principalaxes $sel]
set A [orient $sel [lindex $I 2] {0 0 1}]
$sel move $A
set I [draw principalaxes $sel]
set A [orient $sel [lindex $I 1] {0 1 0}]
$sel move $A
set I [draw principalaxes $sel]

이것을 실행하면 분자의 가장 긴 축이 z-축에 정렬한 것을 볼 수 있다.

Upload:helix.jpg

그런 후 이 분자를 가지고 x와 y축 움직임은 고정시키고 z축으로만 움직이는 MD계산을 하려면 configure file에 다음과 같이 추가한다.

constraints on
consexp 2
consref PDB_FILE_NAME
conskfile PDB_FILE_NAME
conskcol B
selectConstraints on
selectConstrX on
selectConstrY on

이 때, PDB_FILE_NAME은 계산을 하고자 하는 분자의 PDB 파일 이름을 쓰기도 하고 constraint 정보가 있는 다른 PDB 파일을 만들어서 불러올 수도 있다. 여기에서는 같은 파일을 사용한다. conskcol B 는 B column에 constraint constant정보가 있다는 뜻이므로 B column에 constant를 써 넣는다. (B column은 PDB 파일에서 xyz coordinate정보가 나오고 그 다음 2번째 줄에 해당한다. 사실이 컬럼은 temperature B-factor를 쓰는 곳이다.) 여기에서

selectConstrX on
selectConstrY on

대신에

selectConstrZ on

를 쓰면 z-축으로 constraint를 주고 x-y평면으로만 움직이는 MD계산을 할 수 있다.

8 Solvation[ | ]

단백질의 MD나 Monte Carlo계산을 할 때 gas phase에서 하기도 하지만 최근의 경향은 단백질이 물 안에서 보이는 성질을 알고자 하는 것이기 때문에, 물 box를 만드는것은 필수적이다. NAMD에서는 TIP3P model을 지원한다. 여기서, 혼동하지 말아야 할 사실 하나는, 계산은 NAMD에서 하지만 물박스는 VMD에서 만들어야 한다는 것이다.

실제, 물 속에 단백질을 넣는 작업을 실행해 보자. 일단 VMD를 켜고 물 속에 넣기를 원하는 분자를 열어 보자. 물론 이때 psf와 pdb 파일을 모두 가지고 있어야 한다. 만일 pdb 파일만 가지고 있으면 이 NAMDandVMD 시리즈의 맨 처음 글로 돌아가서 psf 파일 만드는 법을 읽어 보고 psf 파일부터 만들어야 한다. 그런 다음 VMD의 command 창으로 가서 다음과 같이 입력한다 :

package require solvate

그런 다음 solvate 라고 치면

 Usage: solvate <psffile> <pdbfile> <option1?> <option2?>...

 Options:
    -o <output prefix> (data will be written to output.psf/output.pdb)
    -s <segid prefix> (should be either one or two letters; default WT)
    -b <boundary> (minimum distance between water and solute, default 2.4)
    -minmax <minmax as returned by measure minmax>
    -t <pad in all directions> (override with any of the following)
    -x <pad negative x>
    -y <pad negative y>
    -z <pad negative z>
    +x <pad positive x>
    +y <pad positive y>
    +z <pad positive z>

과 같이 solvate 명령에 대한 간단한 manual이 나타난다. 이 중에서 -t 옵션을 사용해서 물박스를 만들어보자.

solvate test.psf test.pdb -t 10

라고 입력하면 test.psf와 test.pdb로 만들어진 분자의 주위에 물분자를 뿌리면서 물 박스를 만들게 되는데, 위와 같이 하면 +x, +y, +z, -x, -y, -z의 모든 방향으로 분자와 물박스의 맨 바깥쪽 면과의 거리가 10 옹스트롬이 되도록 물박스를 만든다.

Upload:solvation.jpg

이런 물박스를 갖고 MD 계산을 하려면 periodic boundary condition을 주어야 하고, 일반적으로는 Ewald summation 방법으로 charge를 계산하는 것이 대부분이다. 그렇게 하기 위해서 configure file에 다음과 같은 부분이 필요하다 :

#PBC
cellBasisVector1 49.99 0 0
cellBasisVector2 0 49.99 0
cellBasisVector3 0 0 49.90
cellOrigin 0.810 0.137 -0.899
wrapWater on
wrapAll on
wrapNearest on
#PME
PME on
PMEGridSizeX 32
PMEGridSizeY 32
PMEGridSizeZ 32

3개의 cellBasisVector 명령어는 x, y, z 각 방향으로의 물박스 크기이고 cellOrigin은 물박스의 중심을 의미한다. 이것을 측정하고 싶을 때는 일단 물박스를 만든 후 아래와 같이 하면 구체적인 수치를 알아낼 수 있다 :

set wt1 [atomselect top "segid WT1"]
measure minmax $wt1
measure center $wt1

3개의 PMEGridSize명령은 Ewald summation 계산을 할 때 사용하는 포인트의 숫자를 나타낸다. 이것을 크게 할수록 계산량이 많아지므로 당연히 속도가 느려진다.

9 Input File[ | ]

전형적인 input file은 아래와 같다. 아래에서는 18mer_solvate.psf과 18mer_solvate.pdb의 두 개의 파일을 사용하여 계산을 시작하고, 전체 계산 시간은 1 ns이며 1fs time step을 사용한다. 또한, PME 방법과 periodic boundary condition을 사용한다. 그리고, Langevin dynamics와 Langevin piston을 사용하여 NpT ensemble을 만든다. (더 자세한 설명은 NAMD user manual을 참조)

# NAMD configuration file
# made by One-Sun Lee
# Department of Chemistry, University of Pennsylvania
# 3/7/2003
structure 18mer_solvate.psf
coordinates 18mer_solvate.pdb

set temperature 400.0
temperature $temperature
seed 382730

outputName anneal400k

restartfreq 100
dcdfreq 100
xstFreq 100

outputEnergies 100
outputTiming 100
binaryoutput no
# force field
paratypecharmm on
parameters /home/onesun/toppar/nonbio.prm

exclude scaled1-4
1-4scaling 1.0
switching on
switchdist 8
cutoff 10
pairlistdist 12
stepspercycle 5
#integrator
timestep 1.0
#PBC
cellBasisVector1 39.972 0 0
cellBasisVector2 0 39.940 0
cellBasisVector3 0 0 39.935
cellOrigin 0.09 0.12 -0.14
wrapWater on
wrapAll on
wrapNearest off
#PME
PME on
PMEGridSizeX 40
PMEGridSizeY 40
PMEGridSizeZ 40
#Langevin
langevin on
langevinDamping 5.0
langevinTemp $temperature
langevinHydrogen no
#Nose-Hoover
useGroupPressure yes
useFlexibleCell no
useConstantArea no

LangevinPiston on
LangevinPistonTarget 1.01325
LangevinPistonPeriod 100.0
LangevinPistonDecay 50.0
LangevinPistonTemp $temperature
#minimization on
minimize 10
reinitvels $temperature
run 100000

참조를 위해서 NAMD Tutorial을 방문해 보라. 이 tutorial은 NAMD에서 공식적으로 제공하는 NAMD manual보다 설명이 약간 더 자세하고 예제가 많다.

10 그래핀 옥사이드 만들기[ | ]

PATCH 만들기.

PRES COC         0.00 
GROUP                  
ATOM 1C    CC      0.24 
GROUP
ATOM OS   OS     -0.48 
GROUP
ATOM 2C    CC      0.24 
BOND 1C OS  OS 2C
IC -C  +C   1C     OS    1.4000 120.0000   60.0000 100.0000  1.2000
IC -C  +C   2C     OS    1.4000 100.0000   60.0000 120.0000  1.2000

PRES COC2        0.00 
GROUP                  
ATOM 1C    CC     0.24
GROUP
ATOM OS   OS     -0.48 
GROUP
ATOM 2C    CC     0.24
BOND 1C OS  OS 2C       
IC -C  +C   1C     OS    1.4000 120.0000  -60.0000 100.0000  1.2000
IC -C  +C   2C     OS    1.4000 100.0000  -60.0000 120.0000  1.2000

위 패치를 topology 파일에 써 넣고 VMD 창에서 source topology_name.inp 이라고 한뒤 아래 스크립트 실행.

readpsf gra.psf
coordpdb gra.pdb
set ttt {}
set limit 166

mol load psf gra.psf pdb gra.pdb
set all [atomselect top "all"]
set num [$all num]

proc myrand {num} {
	return [expr int($num*rand())]
}

while {[llength $ttt] < $limit} {
	set atom1 [myrand $num]
	puts "atom1 resid = $atom1"
	if {[lsearch $ttt $atom1]<0} {
		lappend ttt $atom1

		set atom2all [atomselect top "(within 2 of resid $atom1) and (not resid $atom1)"]
		set rand2 [myrand [$atom2all num]]
		set atom2 [lindex [$atom2all get resid] $rand2]
		if {[lsearch $ttt $atom2]<0} {
			lappend ttt $atom2
			if {[expr $atom2%2]==0} {
				patch COC CCC:$atom1 CCC:$atom2
			} else {
				patch COC2 CCC:$atom1 CCC:$atom2
			}
		}
	}
}

guesscoord

writepsf ttt8.psf
writepdb ttt8.pdb

mol delete all
mol load psf ttt8.psf pdb ttt8.pdb

11 GROMACS note[ | ]

기존 분자에 다른 분자를 추가한 시스템을 만드는 경우 (1 mol1 + 2 mol2)

1. gmx insert-molecule -f mol1.gro -ci mol2.gro -o complex.gro -nmol 2 -box 24 24 24

2. gmx solvate -cp complex.gro -cs water.gro -o waterbox.gro -box 24 24 24 -radius 0.3

3. gmx solvate -cp waterbox.gro -cs wf.gro -o wfbox.gro -maxsol #n (#n은 WF의 숫자. 대충 W의 10%로 만든다)

(WF를 넣을때는 water.gro에서 카피해서 wf.gro를 만든다. 이때 column에 주의. W를 모두 WF로 바꿀때 (W --> WF) 원래의 W자리에 F가 오게 만들어야 한다. 즉 원래 W와 나중 F가 같은 column)

4. 이렇게 만든 후 wfbox.gro에서 W와 WF의 숫자를 세어서 top file에 써 넣는다.

5. gmx grompp -f em.mdp -c wfbox.gro -p test.top -maxwarn 10

6. gmx mdrun -v -c minimized.gro

7.gmx grompp -f md_pme.mdp -c minimized.gro -p test.top -maxwarn 10

8. gmx mdrun -v -c

12 같이 보기[ | ]

NAMD, VMD