Amber in practice (in Chinese)

http://www.zhuaxia.com/item/211097608:

第一步：生成小分子模板

蛋白质中各氨基酸残基的力参数是预先存在的，但是很多模拟过程会涉及配体分子，这些有机小分子有很高的多样性，他们的力参数和静电信息不可能预存在库文件中，需要根据需要自己计算生成模板。amber中的antechamber 程序就是生成小分子模板的。

生成模板要进行量子化学计算，这一步可以由antechamber中附带的mopac完成，也可以由gaussian完成，这里介绍用gaussian计算的过程。

建议在计算前用sybyl软件将小分子预先优化，不要用gaussian优化，大基组从头计算进行几何优化花费的计算时间太长。gaussian计算的输入文件可以用antechamber程序直接生成，生成后去掉其中关于几何优化的参数即可

将小分子优化后的结构存储为mol2各式，上传到工作目录，用antechamber程序生成gaussian输入文件，命令如下：

antechamber -i 49.mol2 -fi mol2 -o 49.in -fo gzmat

这样可以生成49.in文件，下载到windows环境，运行gaussian计算这个文件，如果发现计算时间过长或者内存不足计算中断，可以修改文件选择小一些的基组。

获得输出文件49.out之后将它上传到工作目录，再用antechamber生成模板，命令如下：

antechamber -i 49.out -fi gout -o 49mod.mol2 -fo mol2 -c resp

运行之后就会生成一个新的mol2文件，如果用看图软件打开这个文件会发现，原子的颜色很怪异，这是因为mol2的原子名称不是标准的原子名称，看图软件无法识别。下面一步是检查参数，因为可能会有一些特殊的参数在gaff中不存在需要程序注入，命令如下：

parmchk -i 49mod.mol2 -f mol2 -o 49mod

这样那些特殊的参数就存在49mod这个文件中了

第二步：处理蛋白质文件

amber自带的leap程序是处理蛋白质文件的，他可以读入PDB各式的蛋白质文件，根据已有的力场模板为蛋白质赋予键参数和静电参数。

PDB格式的文件有时会带有氢原子和孤对电子的信息，但是在这种格式下氢原子和孤对电子的命名不是标准命名，力场模板无法识别这种不标准的命名，因此需要将两者的信息删除

ATOM 12 1H ARG A 82 12.412 8.891 34.128 1.00 0.00 H

在PDB各式里面，氢原子的信息会在第13或者14列出现H字符，可以应用grep命令删除在13或者14列出现H的行

命令如下：

grep -v '^.............H' 1t4j.pdb > x
grep -v '^............H' x > 1t4j_noh.pdb

除了删除氢和孤对电子，还应该把文件中的结晶水、乙酸等分子删除，这些分子的信息常常集中在文件的尾部，可以直接删除。

处理过之后的蛋白质文件，只包括各氨基酸残基和小分子配体的重原子信息，模拟需要的氢原子和水分子将在leap中添加

接下来需要进一步整理蛋白质文件，主要关注氨基酸的不同存在形态和小分子的原子名称。

半胱氨酸有两种存在形态，有的是两个半胱氨酸通过二硫键相连，有的则是自由的，两种半胱氨酸在力场文件中是用不同的unit来表示的，这相当于是两个完全不同的氨基酸，需要手动更改蛋白质文件中半胱氨酸的名字，桥连的要用CYX，自由的用CYS

组氨酸有若干种质子态，和半胱氨酸一样，也需要查阅文献确定这些质子态，并更改残基名称

最后需要修改的是配体分子的原子名，这是工作量最大的事情，仔细观察可以发现，一般PDB文件中配体的各个原子名称，和我们上面通过antechamber 生成的49mod.mol2中原子名称并不一致，这会造成leap无法识别这些原子，无法为这些原子赋予力参数和静电参数，因此需要手动更改蛋白质文件中配体分子的原子名称。

进行这一步可以同时用看图软件打开49mod.mol2和蛋白质文件，隐藏蛋白质文件中除配体分子以外的所有结构，旋转两个文件，让他们姿态相近，以方便观察，并且在各自均标识原子名称，然后用文字编辑软件打开蛋白质文件，核对看图软件中两个分子对应的原子名称，手动逐一修改。

修改原子名称需要极大的耐心和细心，如果发生错误下一步的操作会无法继续。我现在想到的也只有这个笨办法，如果谁还有别的好法子，欢迎告诉我。

现在文件的准备工作都已经完成，该开始正式的模拟了
第三步：生成拓扑文件和坐标文件

用amber进行分子动力学模拟需要坐标和拓扑文件，坐标文件记录了各个质点所座落的坐标，拓扑文件记录了整个体系各质点之间的链接状况、力参数电荷等信息。这两个文件是由leap 程序生成的

amber中有两个leap程序，一个是纯文字界面的tleap，一个是带有图形界面的Xleap。但是amber的图形界面做得很差，用远程登录使用图形界面不仅麻烦而且慢，所以我比较偏爱使用tleap，两个leap的命令是完全一样的，其实用哪一个都无所谓。

启动tleap，在shell里输入命令tleap，leap就启动了，shell里会显示
-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/prep to search path.
-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/lib to search path.
-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/parm to search path.
-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/cmd to search path.

Welcome to LEaP!
(no leaprc in search path)
>
这个>是leap的提示符

下面要调入库文件。amber是模拟生物分子的好手，主要就是依靠专门为蛋白质多糖核酸量身订做的amber力场，力场的所有参数都存储在库文件里，所以打开leap第一件事便是调入库文件。

amber提供了很多种库，这里我们只用到两个库，gaff和02库，输入命令：
>source leaprc.gaff
>source leaprc.ff02
之后两个库就调入进来了
这时可以用list命令看看库里都有什么：
> list
ACE ALA ARG ASH ASN ASP CALA CARG
CASN CASP CCYS CCYX CGLN CGLU CGLY CHCL3BOX
CHID CHIE CHIP CHIS CILE CIO CLEU CLYS
CMET CPHE CPRO CSER CTHR CTRP CTYR CVAL
CYM CYS CYX Cl- Cs+ DA DA3 DA5
DAN DC DC3 DC4 DC5 DCN DG DG3
DG5 DGN DT DT3 DT5 DTN GLH GLN
GLU GLY HID HIE HIP HIS HOH IB
ILE K+ LEU LYN LYS Li+ MEOHBOX MET
MG2 NALA NARG NASN NASP NCYS NCYX NGLN
NGLU NGLY NHID NHIE NHIP NHIS NILE NLEU
NLYS NMABOX NME NMET NPHE NPRO NSER NTHR
NTRP NTYR NVAL Na+ PHE PL3 POL3BOX PRO
RA RA3 RA5 RAN RC RC3 RC5 RCN
RG RG3 RG5 RGN RU RU3 RU5 RUN
Rb+ SER SPC SPCBOX T4E THR TIP3PBOX TIP4PBOX
TIP4PEWBOXTP3 TP4 TP5 TRP TYR VAL WAT
gaff parm99
>
这里面罗列的就是库里面的unit，包括20种氨基酸、糖以及核酸还有一些常见离子的参数

下面一步是调入配体分子的模板，首先补全参数，输入命令：
>loadamberparams 49mod
然后读入模板文件,输入命令：
>MOL = loadmol2 49mod.mol2
其中MOL是unit的名字，要保证这个名字和pdb文件中配体的残基名完全一致，否则系统仍然无法识别pdb文件中的小分子

现在再输入list命令会发现库里面多了一个unit:
>list
ACE ALA ARG ASH ASN ASP CALA CARG
CASN CASP CCYS CCYX CGLN CGLU CGLY CHCL3BOX
CHID CHIE CHIP CHIS CILE CIO CLEU CLYS
CMET CPHE CPRO CSER CTHR CTRP CTYR CVAL
CYM CYS CYX Cl- Cs+ DA DA3 DA5
DAN DC DC3 DC4 DC5 DCN DG DG3
DG5 DGN DT DT3 DT5 DTN GLH GLN
GLU GLY HID HIE HIP HIS HOH IB
ILE K+ LEU LYN LYS Li+ MEOHBOX MET
MG2 MOL NALA NARG NASN NASP NCYS NCYX
NGLN NGLU NGLY NHID NHIE NHIP NHIS NILE
NLEU NLYS NMABOX NME NMET NPHE NPRO NSER
NTHR NTRP NTYR NVAL Na+ PHE PL3 POL3BOX
PRO RA RA3 RA5 RAN RC RC3 RC5
RCN RG RG3 RG5 RGN RU RU3 RU5
RUN Rb+ SER SPC SPCBOX T4E THR TIP3PBOX
TIP4PBOX TIP4PEWBOXTP3 TP4 TP5 TRP TYR VAL
WAT gaff parm99
>
那个就是配体分子的模板

下面读入pdb文件，输入命令：
>comp = loadpdb 1t4j_noh.pdb
如果输入这个命令之后，屏幕上出现了大量的创建unit或者atom的信息，如下所示，则说明上面一步的pdb文件处理没有做好，还需要重新处理，通常这种情况都发生在配体分子上面，有时则是因为蛋白质中存在特殊残基。
Creating new UNIT for residue: FRJ sequence: 1
Created a new atom named: O36 within residue: .R
Created a new atom named: S33 within residue: .R
Created a new atom named: O35 within residue: .R
Created a new atom named: N34 within residue: .R
如果屏幕仅仅显示添加原子，这说明输入的PDB文件中缺失了部分重原子，leap根据模板自动补全了这些缺失的原子，这种情况不会影响进一步的计算
Added missing heavy atom: .R.A
Added missing heavy atom: .R.A
Added missing heavy atom: .R.A
Added missing heavy atom: .R.A

根据体系的具体情况，还可能要将成对的CYX残基中的二硫键相连，有时候还会链接其他原子，比如将糖基链接在氨基酸残基上，用bond命令可以完成，命令用法如下：

>bond comp.35.SG comp.179.SG

其中comp是刚才读入的分子名称，35和179是残基序号，SG是CYX残基模板中硫原子的名称，用comp.35.SG这样的语法就可以定位一个原子

如果希望进行考虑溶剂效应的动力学模拟，可能还需要为体系加上水，加水有很多种命令，这里只列举一个：

>solvatebox comp TIP3PBOX 10.0

solvatebox命令是说要加上一个方形的周期水箱，comp指要加水的分子，TIP3PBOX是选择的水模板名称，10.0是水箱子的半径

有的体系总电荷不为0，为了模拟稳定，需要加入抗衡离子，系统会自动计算体系的静电场分布，在合适的位置加上离子，命令如下：

>addions comp Na+ 0

意思是用钠离子把体系总电荷补平，还可以选择其他库里面有的离子。

完成这一步之后就可以输出拓扑文件和坐标文件了，但是为了方便起见，在运行最后一步之前要先把leap里加工好的分子单独存成一个库文件，以后可以直接调用这个库文件，免得重复上面的操作：

>saveoff comp 1taj.off

这样就生成了一个off文件在那里，下面输出拓扑文件和坐标文件

>saveamberparm comp 1t4j.prmtop 1t4j.inpcrd
Checking Unit.

Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
total 1 improper torsion applied
Building H-Bond parameters.
Not Marking per-residue atom chain types.
Marking per-residue atom chain types.
(Residues lacking connect0/connect1 -
these don't have chain types marked:

res total affected

CMET 1
)
(no restraints)

>quit

现在准备好拓扑文件和坐标文件，接下来就可以开始能量优化和动力学模拟了。如果愿意的话还可以用ambpdb这个命令生成一个pdb文件，直观地看一看生成了什么东西，命令如下：
[snowyowls@localhost actualamber]$ ambpdb -p 1t4j.prmtop <1t4j.inpcrd> kankan.pdb
| New format PARM file being parsed.
| Version = 1.000 Date = 09/08/06 Time = 16:36:09
[snowyowls@localhost actualamber]$

可以下载之后用看图软件欣赏，如果加了溶剂盒子的话，看的时候可要小心，会恨吓人的样子。

下面就要开始真正的模拟了

转好文章(amber 动力学)

amber 动力学过程：

1、选择力场 tleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99

2、导入晶体结构 model=loadpdb "sbp_lin.pdb"

保存crd和top文件 saveamberparm model polyAT_vac.top polyAT_vac.crd

此时注意电荷是否平衡：

如果缺正电荷 addions model Na+ 0 负离子就用Cl-

选择水箱 solvateoct model TIP3PBOX 8.0

3、保存crd和top文件 saveamberparm model model_wat.top model_wat.crd

4、退出tleap quit

5、保存新的pdb ambpdb -p model_wat.top <> model2.pdb

6、溶剂环境能量最优化。这一步保持溶质（蛋白）不变，去除溶剂中能量不正常的范德华相互作用。

该步骤的配置文件min1.in如下：

---------------------------------------------------------------------

oxytocin: initial minimisation solvent + ions ##说明信息######

&cntrl ##模拟参数起始

imin = 1, ##任务是优化，0 是分子动力学

cut = 10 ##非键相互作用的截断值为10 挨

ntb = 1, ##周期边界条件 0 不采用；1 定容 ；2 定压

maxcyc = 4000, ##优化步数

ntr = 1, ##优化时需要一些约束原子 -ref

ncyc = 2000, ##前2000最陡下降，后面步骤共轭梯度

/

Hold the protein fixed ##约束说明

500.0 ##作用在肽键上的力 kcal/mol

RES 1 9 ##限制的残基序号 同restrain=’:1-9’

END

END

------------------------------------------------------------------------------

任务命令：如果

sander -O -i min1.in -p model_wat.top -c model_wat.crd -o min1.out -r min1.rst –ref model_wat.crd &

7、对蛋白进行优化，min2.in文件将min1.in中的限制原子修改，限制水的位置。

也可以考虑利用restrainmask=’:1-9@CA,N,C’约束蛋白主链上的原子。

sander -O -i min2.in -p model_wat.top -c min1.rst -o min2.out -r min2.rst&

8、整体的优化,去掉限制条件

sander -O -i min3.in -p model_wat.top -c min2.rst -o min3.out -r min3.rst &

9、有限制的分子动力学

第一步分子动力学保持蛋白分子位置不变，但是不是完全固定每个原子，同时缓解蛋白分子周围的水分子，是溶剂环境能量优化。在这个步骤中，我们将主要目的是对特定的原子使用作用力使其能量优化。

Eq1.in 如下：

fix protein ,relax H2O

&cntrl

nstlim=25000, dt=0.002, ntx=1, irest=0, ntpr=500, ntwr=500,ntwx=500,

tempi=0.0,temp0=300,ntt=3, gamma_ln=1.0,

ntb=1, ntp=0,

nrespa=1,

cut = 10,

ntc=2,ntf=2,

NTR=1,

/

fix protein and HEM

10

RES 1 284

END

END

-----------------------------

nstlim = #：＃表示计算的步数。dt = 0.002：表示步长，单位为ps，0.002表示2fs。ntx=1 irest=0 默认 ntb = 1：表示分子动力学过程保持体积固定。

imin = 0：表示模拟过程为分子动力学，不是能量最优化。

temp0 = 300：表示最后系统到达并保持的温度，单位为K。

tempi = 100：系统开始时的温度。 ntc=2,ntf=2 忽略氢键

gamma_ln = 1：表示当ntt＝3时的碰撞频率，单位为ps-1（请参考AMBER手册）

ntt = 3：温度转变控制，3表示使用兰格氏动力学。

sander -O -i eq1.in -p model_wat.top -c min3.rst -o eq1.out -r eq1.rst -ref min3.rst -x eq1.mdcrd

11整系统分子动力学模拟: eq2

-------------------------------------------------

f2:500ps MD

&cntrl

imin = 0, irest=1, ntx=5,

ntb=2, pres0 = 1.0, ntp=1,

taup = 2.0, ntc=2, ntf=2,

cut = 10, ntr = 0,

ntt = 3, gamma_ln = 1.0,

tempi = 300.0 , temp0 = 300.0

nstlim=500000, dt=0.002,

ntpr=500, ntwr=500, ntwx=500

/

----------------------------------------------------------

ntb=2：表示分子动力学过程的压力常数。Pres0＝1：引用1个单位的压强。

ntp＝1：表示系统动力学过程各向同性。taup = 2.0：压力缓解时间，单位为ps。

使用以下命令进行MD：

sander -O -i eq2.in -p model_wat.top -c eq1.rst -o eq2.out -r eq2.rst -x eq2.mdcrd -ref eq1.rst &

结果处理：

在动力学过程中，将会产生 .rst 和.mdrcd 文件（需要的），由于crd文件很多，可以将其压缩成.gz 文件：gzip filename，将产生一个filename.gz 文件

做出已跑动力学的rmd图，判断是否平衡：

ptraj xxx.top <>

xxx.in 内容：

trajin eq2.mdcrd.gz

trajin eq3.mdcrd.gz

trajin eq4.mdcrd.gz

trajin eq5.mdcrd.gz

rms first mass out xin.rms.dat1 :1-284@CA,C,N time 0.1

#产生一个xin.rms.dat1的文件，整体1-284骨架上的C与N原子

rms first mass out xin.rms.dat2 :88-91,172,201-204,230@CA,C,N time 0.1

#产生一个xin.rms.dat2的文件，保守残基骨架上的C与N原子

----------------------------------------------------------------------------

利用xmgrace 作图：

xmgrace xin.rms.dat1 xin.rms.dat2

如果需要取随即的点的构型：

ptraj xxx.top <>

xxx.in内容：

trajin eq7.mdcrd.gz 117 117 #eq7中的117ps 的构型

strip :WAT #冒号前有空格，后没有，注意wat与top的水的大小写一致

strip :Na+

trajout eq7.pdb pdb nobox 产生一个eq7.pdb.117的文件

-------------------------------------------

traj AR.top <<>

> trajin AR.md7.crd 300 500

> center :1-250

> image center familiar

> rms first mass out rms.dat :1-250

> average average.rst restrt

> average average.pdb pdb

> EOF

对于有些小分子利用sybyl构建的一般缺少原子参数，首先将小分子单独存成mol2，用amber的antechamber 做出参数。

1、利用gaussian

* antechamber -i 49.mol2 -fi mol2 -o 49.in -fo gzmat

这样可以生成49.in文件，下载到windows环境，运行gaussian计算这个文件，如果发现计算时间过长或者内存不足计算中断，可以修改文件选择小一些的基组。

you have Gassian output file.out (don’t forget to put some keyword in the Gassian input for the special format used by amber)

* -antechamber –i file.out –fi gout –o filr.prep –fo prepi –c resp (this step prepare the “file.prep” for your drug -c bcc all right )

* -parmchk –i file.prep –f prepi –o file.frcmod )this step prepare the file “file.frcmod”for the drug)

2、 直接用mol2 转

也可以直接把sybyl 的mol2 做参数

* -antechamber –i file.mol2 –fi mol2 –o filr.prep –fo prepi –c resp / -c bcc

* -parmchk –i file.prep –f prepi –o file.frcmod

now you have 2 files: file.prep and file.frcmod for the drug .

注意修改file.prep 中文件名。要与pdb中的配体一致。（vi prep。假设是abc）

对应pdb与prep文件中的原子名称,，原子顺序要一致。

prep 在xleap中打开

using tleap for the protein or drug or/and complex

-xleap –s –f leaprc.ff99

-mod = loadamberparams file.frcmod

-loadamberprep file.prep

edit abc

pdb 可以用sybyl打开。用vi 修改pdb

-source leaprc.gaff

-RL = loadpdb file.pdb

-Solvateoct RL TIP3PBOX 10

-Addions RL Cl- 0 (0:zero for neutral charge,Cl- can be replaced by Na+)

-Saveamberparm RL fileWT.top fileWT.crd (save topology and coord. With water)

Amber网上教程地址[zz]

http://www.rosswalker.co.uk/tutorials/amber_workshop/ http://amber.scripps.edu/tutorials/ amber8/9 http://amber.scripps.edu/dbase.html amber ff in gmx http://folding.stanford.edu/ffamber/ other http://www.chpc.utah.edu/%7Echeatham/software.html AMBER 上好用低SMD计算 各位好以文字附上AMBER低输入档( umbrella.RST # Change distance between atoms 1771 and 3787 from 63 A to 253 A &rst iat=1771,3787,r1=10., r2=63.,r3=63.,r4=210.,rk2=10.,rk3=10., / smd.in Sample pulling input &cntrl nstlim=1000000, cut=99.0, igb=1, saltcon=0.1, ntpr=10, ntwr=10, ntt=3, gamma_ln=5.0, ntx=5, irest=1, ntwx=10, ig = 256251, ntc=2, ntf=2, tol=0.000001, dt=0.002, ntb=0, tempi=310., temp0=310., jar=1, / &wt type='DUMPFREQ', istep1=1, / &wt type='END', / DISANG=dist.RST DUMPAVE=dist_vs_t LISTIN=POUT LISTOUT=POUT amber中生成小分子模板 amber, 小分子, 模板 第一步：生成小分子模板 蛋白质中各氨基酸残基的力参数是预先存在的，但是很多模拟过程会涉及配体分子，这些有机小分子有很高的多样性，他们的力参数和静电信息不可能预存在库文件中，需要根据需要自己计算生成模板。amber中的antechamber 程序就是生成小分子模板的。生成模板要进行量子化学计算，这一步可以由antechamber中附带的mopac完成，也可以由gaussian完成，这里介绍用gaussian计算的过程。 建议在计算前用sybyl软件将小分子预先优化，不要用gaussian优化，大基组从头计算进行几何优化花费的计算时间太长。gaussian计算的输入文件可以用antechamber程序直接生成，生成后去掉其中关于几何优化的参数即可 将小分子优化后的结构存储为mol2各式，上传到工作目录，用antechamber程序生成gaussian输入文件，命令如下： antechamber -i 49.mol2 -fi mol2 -o 49.in -fo gzmat 这样可以生成49.in文件，下载到windows环境，运行gaussian计算这个文件，如果发现计算时间过长或者内存不足计算中断， 可以修改文件选择小一些的基组。 获得输出文件49.out之后将它上传到工作目录，再用antechamber生成模板，命令如下： antechamber -i 49.out -fi gout -o 49mod.mol2 -fo mol2 -c resp 运行之后就会生成一个新的mol2文件，如果用看图软件打开这个文件会发现，原子的颜色很怪异，这是因为mol2的原子名称 不是标准的原子名称，看图软件无法识别。下面一步是检查参数，因为可能会有一些特殊的参数在gaff中不存在需要程序注入， 命令如下： parmchk -i 49mod.mol2 -f mol2 -o 49mod 这样那些特殊的参数就存在49mod这个文件中了 第二步：处理蛋白质文件 amber自带的leap程序是处理蛋白质文件的，他可以读入PDB各式的蛋白质文件，根据已有的力场模板为蛋白质赋予键参数和静电参数。 PDB格式的文件有时会带有氢原子和孤对电子的信息，但是在这种格式下氢原子和孤对电子的命名不是标准命名，力场模板无法识别这 种不标准的命名，因此需要将两者的信息删除 ATOM 12 1H ARG A 82 12.412 8.891 34.128 1.00 0.00 H 在PDB各式里面，氢原子的信息会在第13或者14列出现H字符，可以应用grep命令删除在13或者14列出现H的行 命令如下： grep -v '^.............H' 1t4j.pdb > x grep -v '^............H' x > 1t4j_noh.pdb 除了删除氢和孤对电子，还应该把文件中的结晶水、乙酸等分子删除，这些分子的信息常常集中在文件的尾部，可以直接删除。 处理过之后的蛋白质文件，只包括各氨基酸残基和小分子配体的重原子信息，模拟需要的氢原子和水分子将在leap中添加 接下来需要进一步整理蛋白质文件，主要关注氨基酸的不同存在形态和小分子的原子名称。 半胱氨酸有两种存在形态，有的是两个半胱氨酸通过二硫键相连，有的则是自由的，两种半胱氨酸在力场文件中是用不同的unit来表示的， 这相当于是两个完全不同的氨基酸，需要手动更改蛋白质文件中半胱氨酸的名字，桥连的要用CYX，自由的用CYS 组氨酸有若干种质子态，和半胱氨酸一样，也需要查阅文献确定这些质子态，并更改残基名称 最后需要修改的是配体分子的原子名，这是工作量最大的事情，仔细观察可以发现，一般PDB文件中配体的各个原子名称，和我们上面通过 antechamber 生成的49mod.mol2中原子名称并不一致，这会造成leap无法识别这些原子，无法为这些原子赋予力参数和静电参数， 因此需要手动更改蛋白质文件中配体分子的原子名称。 进行这一步可以同时用看图软件打开49mod.mol2和蛋白质文件，隐藏蛋白质文件中除配体分子以外的所有结构，旋转两个文件， 让他们姿态相近，以方便观察，并且在各自均标识原子名称，然后用文字编辑软件打开蛋白质文件，核对看图软件中两个分子对 应的原子名称，手动逐一修改。 修改原子名称需要极大的耐心和细心，如果发生错误下一步的操作会无法继续。我现在想到的也只有这个笨办法，如果谁还有别的好法子 ，欢迎告诉我。 现在文件的准备工作都已经完成，该开始正式的模拟了 第三步：生成拓扑文件和坐标文件 用amber进行分子动力学模拟需要坐标和拓扑文件，坐标文件记录了各个质点所座落的坐标，拓扑文件记录了整个体系各质点之间的 链接状况、力参数电荷等信息。这两个文件是由leap 程序生成的 amber中有两个leap程序，一个是纯文字界面的tleap，一个是带有图形界面的Xleap。但是amber的图形界面做得很差，用远程登录使用 图形界面不仅麻烦而且慢，所以我比较偏爱使用tleap，两个leap的命令是完全一样的，其实用哪一个都无所谓。 启动tleap，在shell里输入命令tleap，leap就启动了，shell里会显示 -I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/prep to search path. -I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/lib to search path. -I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/parm to search path. -I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/cmd to search path. Welcome to LEaP! (no leaprc in search path) > 这个>是leap的提示符 下面要调入库文件。amber是模拟生物分子的好手，主要就是依靠专门为蛋白质多糖核酸量身订做的amber力场，力场的所有参数都存 储在库文件里，所以打开leap第一件事便是调入库文件。 amber提供了很多种库，这里我们只用到两个库，gaff和02库，输入命令： >source leaprc.gaff >source leaprc.ff02 之后两个库就调入进来了 这时可以用list命令看看库里都有什么： 这里面罗列的就是库里面的unit，包括20种氨基酸、糖以及核酸还有一些常见离子的参数 下面一步是调入配体分子的模板，首先补全参数，输入命令： >loadamberparams 49mod 然后读入模板文件,输入命令： >MOL = loadmol2 49mod.mol2 其中MOL是unit的名字，要保证这个名字和pdb文件中配体的残基名完全一致，否则系统仍然无法识别pdb文件中的小分子 现在再输入list命令会发现库里面多了一个unit: 那个就是配体分子的模板 下面读入pdb文件，输入命令： >comp = loadpdb 1t4j_noh.pdb 如果输入这个命令之后，屏幕上出现了大量的创建unit或者atom的信息，如下所示，则说明上面一步的pdb文件处理没有做好， 还需要重新处理，通常这种情况都发生在配体分子上面，有时则是因为蛋白质中存在特殊残基。 Creating new UNIT for residue: FRJ sequence: 1 Created a new atom named: O36 within residue: .R Created a new atom named: S33 within residue: .R Created a new atom named: O35 within residue: .R Created a new atom named: N34 within residue: .R 如果屏幕仅仅显示添加原子，这说明输入的PDB文件中缺失了部分重原子，leap根据模板自动补全了这些缺失的原子，这种情况 不会影响进一步的计算 Added missing heavy atom: .R.A Added missing heavy atom: .R.A Added missing heavy atom: .R.A Added missing heavy atom: .R.A 根据体系的具体情况，还可能要将成对的CYX残基中的二硫键相连，有时候还会链接其他原子，比如将糖基链接在氨基酸残基上 ，用bond命令可以完成，命令用法如下： >bond comp.35.SG comp.179.SG 其中comp是刚才读入的分子名称，35和179是残基序号，SG是CYX残基模板中硫原子的名称，用comp.35.SG这样的语法就可以定 位一个原子 如果希望进行考虑溶剂效应的动力学模拟，可能还需要为体系加上水，加水有很多种命令，这里只列举一个： >solvatebox comp TIP3PBOX 10.0 solvatebox命令是说要加上一个方形的周期水箱，comp指要加水的分子，TIP3PBOX是选择的水模板名称，10.0是水箱子的半径 有的体系总电荷不为0，为了模拟稳定，需要加入抗衡离子，系统会自动计算体系的静电场分布，在合适的位置加上离子，命令如下： >addions comp Na+ 0 意思是用钠离子把体系总电荷补平，还可以选择其他库里面有的离子。 完成这一步之后就可以输出拓扑文件和坐标文件了，但是为了方便起见，在运行最后一步之前要先把leap里加工好的分子单独存 成一个库文件，以后可以直接调用这个库文件，免得重复上面的操作： >saveoff comp 1taj.off 这样就生成了一个off文件在那里，下面输出拓扑文件和坐标文件 >saveamberparm comp 1t4j.prmtop 1t4j.inpcrd Checking Unit. Building topology. Building atom parameters. Building bond parameters. Building angle parameters. Building proper torsion parameters. Building improper torsion parameters. total 1 improper torsion applied Building H-Bond parameters. Not Marking per-residue atom chain types. Marking per-residue atom chain types. (Residues lacking connect0/connect1 - these don't have chain types marked: res total affected CMET 1 ) (no restraints) >quit 现在准备好拓扑文件和坐标文件，接下来就可以开始能量优化和动力学模拟了。如果愿意的话还可以用ambpdb这个命令生成一 个pdb文件，直观地看一看生成了什么东西，命令如下： [snowyowls@localhost actualamber]$ ambpdb -p 1t4j.prmtop <1t4j.inpcrd> kankan.pdb | New format PARM file being parsed. | Version = 1.000 Date = 09/08/06 Time = 16:36:09 [snowyowls@localhost actualamber]$ 可以下载之后用看图软件欣赏，如果加了溶剂盒子的话，看的时候可要小心，会恨吓人的样子。 第四步：能量优化 用amber进行分子动力学模拟需要坐标和拓扑文件，这在上一步已经完成了，分别生成了1t4j.prmtop 和1t4j.inpcrd两个文件 ，下面一步就是动力学模拟之前的能量优化了。 由于我们进行的起始结构来自晶体结构或者同源模建，所以在分子内部存在着一定的张力，能量优化就是在真正的动力学之前， 释放这些张力，如果没有这个步骤，在动力学模拟开始之后，整个分子可能会因此散架。 能量优化由sander模块完成，运行sander至少需要三个输入文件，分别是分子的拓扑文件，坐标文件以及sander的控制文件。 现在分子的拓扑文件和坐标文件已经完成，需要建立输入文件，min_1.in Initial minimisation of our structures &cntrl imin=1, maxcyc=4000, ncyc=2000, cut=10, ntb=1,ntr=1, restraint_wt=0.5 restraintmask=':1-283' / 文件首行是说明，说明这项任务的基本情况； &cntrl与/之间的部分是模拟的参数 其中imin=1表示任务是能量优化，maxcyc=4000表示能量优化共进行4000步，ncyc=2000表示在整个能量优化的4000步中， 前 2000步采用最陡下降法，在2000步之后转换为共轭梯度法，如果模拟的时候不希望进行方法的转换，可以再加入另一 个关键词NTMIN，如果NTMIN =0则全程使用共轭梯度法，NTMIN=2则全程使用最陡下降法，此外还有=3和=4的选项，分别是xmin法和lmod法 ，具体情况可以看手册。 第二行的cut=10表示非键相互作用的截断值，单位是埃， ntb=1表示使用周期边界条件，这个选项要和前面生成的拓扑文件坐标文件相匹配， 如果前面加溶剂时候用的是盒子水，就设置ntb=1，如果加的是层水，那就应该选择ntb=0；ntr=1表示在能量优化的过程中要约束一些原子， 约束的是哪些原子呢？后面有。 第三行和第四行都是关于约束原子的信息，restraint_wt=0.5限定了约束的力常数，在这里约束原子就是把原子用一根弹簧拉在固定的 位置上，一旦原子偏离固定的位置，系统就会给他施加一个回复力，偏离的越远，回复力越大，回复力就是由这个力常数决定的，单位是 Kcal/(mol*A)。 restraintmask=':1-283'表示约束的是从1到283号残基，在这个分子中，1-283号残基是蛋白质上的氨基酸残基，从284 号开始，就都是水了，所以这个命令的意思就是，约束整个蛋白质，自由地优化溶剂分子。因为溶剂分子是前面的tleap自动给加上的， 所以一定要额外多关注一些。 下面运行sander来执行这个优化： [snowyowls@localhost actualamber]$ sander -O -i min_1.in -p 1t4j.prmtop -c 1t4j.i npcrd -ref 1t4j.inpcrd -r 1t4j_min1.rst -o 1t4j_min1.out 命令中，-O表示覆盖所有同名文件，-i min_1.in表示sander的控制文件是min_1.in，-p 1t4j.prmtop表示分子的拓扑文件 ，-c 1t4j.inpcrd表示坐标文件，-ref 1t4j.inpcrd是参考坐标文件，只有在控制文件中出现关键词ntr=1的时候才需要给 sander指定-ref文件，这是约束原子的参考坐标，- ref 1t4j.inpcrd就是说以1t4j.inpcrd中的坐标为准进行约束原子的优化。 以上这四个是输入文件。-r 1t4j_min1.rst表示经过模拟之后新的原子坐标会输出到1t4j_min1.rst文件中，-o 1t4j_min1.out 则表示优化过程中的相关信息都会写入到1t4j_min1.out文件中。 运行起这个命令之后，等拿到 1t4j_min1.rst文件就意味着对溶剂的优化已经差不多了，显然下面还需要对蛋白质本身进行优化， 这个优化还要分两步进行，控制文件分别是min_2.in 和min_3.in min_2.in Initial minimisation of our structures &cntrl imin=1, maxcyc=5000, ncyc=2500, cut=10, ntb=1,ntr=1, restraint_wt=0.5 restraintmask=':1-283@CA,N,C' / 在这里发生变化的是约束原子的范围， ':1-283@CA,N,C'表示1-283号残基中名叫CA，N和C的原子，这些原子实际上是蛋白质主链上 的原子，也就是说这一次的优化是约束了蛋白质主链上的原子之后，对溶剂和侧链原子进行自由优化。 min_3.in Initial minimisation of our structures &cntrl imin=1, maxcyc=10000, ncyc=5000, cut=10, ntb=1, / 在这里不再进行约束原子的优化了，对整个分子进行全原子优化。 三步优化的命令如下： [snowyowls@localhost actualamber]$ sander -O -i min_1.in -p 1t4j.prmtop -c 1t4j.inpcrd -ref 1t4j.inpcrd -r 1t4j_min1.rst -o 1t4j_min1.out [snowyowls@localhost actualamber]$ sander -O -i min_2.in -p 1t4j.prmtop -c 1t4j_min1.rst -ref 1t4j_min1.rst -r 1t4j_min2.rst -o 1t4j_min2.out [snowyowls@localhost actualamber]$ sander -O -i min_3.in -p 1t4j.prmtop -c 1t4j_min2.rst -r 1t4j_heat0.rst -o 1t4j_min3.out 接下来就是真正的分子动力学模拟了……

本文引用地址： http://www.sciencenet.cn/m/user_content.aspx?id=289031

AMBER分子动力学简例[zz]

AMBER分子动力学简例（一） 概述 以下是使用AMBER包的简单教程，希望对开始学习分子动力学的同学有用处。申明一下，以下教程原版来自网上，是最最基本的教程，同时也非常实用，有非常好的借鉴意义。 AMBER分子动力学程序包是加州圣弗兰西斯科大学（University of California San Francisco,UCSF)的Peter A Kollman和其同事编的，程序很全，现在已经发展到版本9.0。AMBER功能涵盖种类非常多的生物分子，包括蛋白、核算以及药物小分子。软件详细情况请浏览http://amber.scripps.edu. 以下是AMBER软件包中四个主要的大程序： -------------------------------------------------------------------------------- Leap：用于准备分子系统坐标和参数文件，有两个程序： xleap：X－windows版本的leap，带GUI图形界面。 tleap：文本界面的Leap。 Antechamber：用于生成少见小分子力学参数文件的。有的时候一些小分子Leap程序不认识， 需要加载其力学参数，这些力学参数文件就要antechamber生成。 Sander：MD数据产生程序，即MD模拟程序，被称做AMBER的大脑程序。 Ptraj：MD模拟轨迹分析程序。 -------------------------------------------------------------------------------- 学习项目 本教程研究的题目是脑下垂体荷尔蒙之一的oxytocin，需要X光衍射晶体结构文件1NPO.PDB。该文件包含了该荷尔蒙和其运载蛋白的复合物，可以从蛋白数据库下载。 PDB文件是不包含氢原子的，Leap程序会自动的加上PDB文件缺少的东西。当第一次使用PDB文件的时候，要十分留意文件包含的信息，所以PDB文件缺少的残疾、侧链或者添加的变异都在这个地方记录。可以用文本阅读程序阅读PDB文件头部的信息，即以REMARKS开头的信息文本行。PDB一个重要的信息是SSBOND记录，该记录说明结构中二硫键的位置，这样的信息在使用Leap程序建立分子结构的时候需要。在本教程中，我们将比较oxytocin在真空和溶液中分子动力学的差异，如果没有二硫键，将影响整个结果。整个过程在Linux系统下完成，如对Linux不熟悉，请翻阅有关Linux书籍。 建立项目目录，并进入该目录： mkdir myproj cd myproj 下载1NPO.PDB文件到该目录下，使用文本阅读器阅读该PDB文件。从在文件的开头部分可以得到该文件是一个二聚物的PDB文件，删除文件中A、C、和D链对应的文本行，剩下的就是oxytocin的晶体机构了，保存为oxyt.pdb文件。 使用Swiss PDB viewer分析以下oxyt.pdb文件，可以得知该pdb文件缺少了一个侧链。如果使用Deep View，它会自动给文件加上侧链。重新保存pdb文件为oxyt.pdb文件。 xleap和tleap程序功能是一样的，都是准备分子结构的坐标文件和拓扑文件。xleap启动一个X界面，比较慢。tleap纯文本，要快一些，我们选择tleap： tleap -s -f leaprc.ff03 其中leaprc.ff03是AMBER的2003力场文件。力场是一个很重要的文件，定义了分子、原子和残疾等的信息，请查阅有关资料。 oxy=loadpdb oxyt.pdb 该命令读入oxyt.pdb文件到oxy变量中。 bond oxy.1.SG oxy.6.SG 连接oxy变量的第一和第六个残疾的SG原子，即是连接二硫键。可以使用check命令检查oxy是否完好，没有特殊情况的话，最后应该是“OK”，表示分子系统状态是好的，可以保存。 check oxy 接下来就是保存分子系统了： saveamberparm oxy oxy_vac.top oxy_vac.crd 其中oxy_vac.top和oxy_vac.crd分别是分子系统真空状态下的拓扑和坐标文件。 接下来要给分子系统添加水环境， solvateoct oxy TIP3PBOX 9.0 该命令让Leap程序使用TIP3PBOX水模型,水环境距离分子为9纳米。也可以用solvatebox命令，但是solvatebox添加的水环境是一个立方体，不是八面体。一般要求水表面到蛋白距离要达到8.5纳米，避免MD过程中蛋白冲出水环境，但是水环境越大计算时间将越长。如果MD过程内含PME，那么水箱长度要大于2倍截矩（cutoff），截矩在MD配置文件中定义。加溶剂之后，要计算系统是否为电中性，使用命令： charge oxy 如果现实不是电中性，要使用addions添加反性电荷，如Cl－或者Na＋等。 最后保存溶剂环境的系统拓扑结构和坐标文件，并退出Leap程序： saveamberparm oxy oxy.top oxy.crd quit 也可以把命令集中到一个，让tleap读取。如将上面的命令集中到以下文件中 "oxy.leaprc" -------------------------------------------------------------------------------- source leaprc.ff03 oxy=loadpdb oxyt.pdb bond oxy.1.SG oxy.6.SG check oxy saveamberparm oxy oxy_vac.top oxy_vac.crd solvateoct oxy saveamberparm oxy oxy.top oxy.crd quit -------------------------------------------------------------------------------- 将以上两横线之间的命令保存为"oxy.leaprc"，然后使用以下命令一步搞定： tleap -s -f oxy.leaprc 到此，分子系统准备已经完成，目录中会多了一个leap.log文件，这是Leap程序的记录文件，它包含了leap程序所做的一切，包括自动的和手动命令。 现在已经有了可以进行分子动力学的基本文件了，再编写sander的控制文件，就可以进行模拟了。这些将在后文中介绍。 AMBER分子动力学简例（二） 分子动力学（1） 真空模式 真空模式分子动力学模拟将使用NVT系宗分两步进行，即系统能量最优化和分子动力学过程。 1、系统能量最优化。我们将使用淬火能量最优化解除系统内原子之间的不正常相互作用，这些原子之间的高能量相互作用如果不消除，可能影响后续的分子动力学过程。因为动力学过程是能量梯度变化的，太高的能量壁垒可能让MD局限在某一个能量局部最小化位置中。 Sander程序是分子动力学模拟程序，它的主要功能是能量最优化，动力学模拟和NMR优化计算。我们必须给Sander一个运行的配置文件，使其按照我们的要求进行计算。能量最优化的配置文件如下： "min_vac.in" -------------------------------------------------------------------------------- oxytocin: initial minimization prior to MD &cntrl imin = 1, maxcyc = 500, ncyc = 250, ntb = 0, igb = 0, cut = 12 / -------------------------------------------------------------------------------- Sander程序的基本命令参数如下： sander –O –i in –o out –p prmtop –c inpcrd –r restrt [-ref refc –x mdcrd –v mdvel –e mden –inf mdinfo] 所以我们的命令可以如下： sander –O –i min_vac.in –o min_vac.out –p oxy_vac.top –c oxy_vac.crd –r oxy_vacmin.rst & 可以使用more命令或者tail命令查看min_vac.out的输出内容，那是能量最优化的记录。 2、分子动力学模拟。 Sander程序的配置文件为： md_vac.in -------------------------------------------------------------------------------- oxytocin MD in-vacuo, 12 angstrom cut off, 250 ps &cntrl imin = 0, ntb = 0, igb = 0, ntpr = 100, ntwx = 500, ntt = 3, gamma_ln = 1.0, tempi = 300.0, temp0 = 300.0, nstlim = 125000, dt = 0.002, cut = 12.0 / -------------------------------------------------------------------------------- 命令为： sander –O –i md_vac.in –o md_vac.out –p oxy_vac.top –c oxy_vacmin.rst –r oxy_vacmd.rst –x oxy_vacmd.mdcrd –ref oxy_vacmin.rst –inf mdvac.info & MD的计算过程一般比较久，真空相对与溶剂中要快以下，250ps的模拟大概在一个主频为2。0GHz的Linux单机上运行5分钟。 以上配置文件中用到很多参数，这些参数这是sander程序参数的一小部分，以下将相应解释。如果要了解sander的其他参数，请阅读AMBER用户指南。 &ctrl和"/"：sander的参数一般要求出现在这两个标识符号之间，参数以及这两个标识符被称做控制模块。 cutoff：以纳米为单位的截矩。即超出截矩范围的非键连接相互作用将不计。 ntr：原子位置能量抑制位，1表示抑制，0表示不抑制。 imin：能量最优化标志位。1表示sander将进行能量最优化，0表示让sander进行分子动力学模拟。 macyc：能量最优化次数。 ncyc：便是经过多少次能量优化以后，能量优化从淬火过程变为梯度变化过程。 ntmin：能量优化方法标志位。0表示前10个能量最优化为淬火过程，然后进行梯度能量优化；1表示ncyc次淬火过程，然后进行梯度能量优化，为默认值；2表示只进行淬火过程。 dx0：表示启动模拟步长。 dxm：最大优化步数。 drms：梯度能量优化标准，默认值为1.0E-4 kcal/mol.A。 更多参数将在后文中解释。 AMBER分子动力学简例（三） 分子动力学（2） 水环境中的分子动力学模拟 溶剂环境中的分子动力学模拟分为以下四步进行： 1、溶剂环境能量最优化。这一步保持溶质（蛋白）不变，去除溶剂中能量不正常的范德华相 互作用。 2、整系统能量最优化。去除整个系统中能量不正常的相互作用。 3、有限制的分子动力学。保持蛋白质不动，溶解溶剂的不同层，同时逐渐将系统温度从0K提升到300K。 4、整系统分子动力学模拟。在一个大气压，300K的环境下整个系统分子动力学模拟。可以得到成果的分子动力学模拟。 ############################# 1、溶剂环境能量最优化。 该步骤的配置文件min1.in如下： -------------------------------------------------------------------------------- oxytocin: initial minimisation solvent + ions &cntrl imin = 1, maxcyc = 1000, ncyc = 500, ntb = 1, ntr = 1, cut = 10 / Hold the protein fixed 500.0 RES 1 9 END END -------------------------------------------------------------------------------- 该过程保持肽链不动，其中500.0单位是kcal/mol，表示作用在肽链上使其不动的力。“RES 1 9”表示肽链残基数目，因为我们学习使用的oxytocin有9个残基。 模拟命令如下： sander –O –i min1.in –o min1.out –p oxy.top –c oxy.crd –r oxy_min1.rst –ref oxy.crd & 2、整系统能量最优化。 配置文件min2.in如下： -------------------------------------------------------------------------------- oxytocin: initial minimisation whole system &cntrl imin = 1, maxcyc = 2500, ncyc = 1000, ntb = 1, ntr = 0, cut = 10 / -------------------------------------------------------------------------------- 命令如下： sander –O –i min2.in –o min2.out –p oxy.top –c oxy_min1.rst –r oxy_min2.rst & 3、有限制的分子动力学。 第 一步分子动力学保持蛋白分子位置不变，但是不是完全固定每个原子，同时缓解蛋白分子周围的水分子，是溶剂环境能量优化。在这个步骤中，我们将主要目的是对 特定的原子使用作用力使其能量优化。我们要优化溶剂环境，至少需要10ps，我们将使用20ps用来优化我们上两步制作的分子系统的周期性边界的溶剂环 境。 命令配置文件md1.in如下： -------------------------------------------------------------------------------- oxytocin: 20ps MD with res on protein &cntrl imin = 0, irest = 0, ntx = 1, ntb = 1, cut = 10, ntr = 1, ntc = 2, ntf = 2, tempi = 0.0, temp0 = 300.0, ntt = 3, gamma_ln = 1.0, nstlim = 10000, dt = 0.002, ntpr = 100, ntwx = 500, ntwr = 1000 / Keep protein fixed with weak restraints 10.0 RES 1 9 END END -------------------------------------------------------------------------------- 上述参数解释如下： ntb = 1：表示分子动力学过程保持体积固定。 imin = 0：表示模拟过程为分子动力学，不是能量最优化。 nstlim = #：＃表示计算的步数。 dt = 0.002：表示步长，单位为ps，0.002表示2fs。 temp0 = 300：表示最后系统到达并保持的温度，单位为K。 tempi = 100：系统开始时的温度。 gamma_ln = 1：表示当ntt＝3时的碰撞频率，单位为ps-1（请参考AMBER手册） ntt = 3：温度转变控制，3表示使用兰格氏动力学。 tautp = 0.1：热浴时间常数，缺省为1.0。小的时间常数可以得到较好的耦联。 vlimit = 20.0：保持分子动力学稳定性速度极限。20.0为缺省值，当动力学模拟中原子速度大于极限值时，程序将其速度降低到极限值以下。 comp = 44.6： 溶剂可压缩单位。 ntc = 2：Shake算法使用标志位。1表示不实用使用，2表示氢键将被计算，3表示所有键都将被计算在内。 tol = #.#####：坐标位置重新设置的几何位置相对容忍度。 我 们将使用一个较小的作用力，10kcal/mol。在分子动力学中，当ntr＝1时，作用力只需要5-10kcal/mol（我们需要引用一个坐标文件做 分子动力学过程的比较，我们需要使用"-ref"参数）。太大的作用力同时使用Shake算法和2fs步长将使整个系统变得不稳定，因为大的作用力使系统 中的原子产生大频率的振动，模拟过程并步需要。 运行命令如下： sander –O –i md1.in –o md1.out –p oxy.top –c oxy_min2.rst –r oxy_md1.rst –x oxy_md1.mdcrd –ref oxy_min2.rst –inf md1.info& 进行MD的运行时间一般较长，可以使用程序的并行版本提交集群计算。在主频位2.0GHz的P4单机上，大概需要一个钟头。可以随时查看md1.in文件的程序输出。 4、整系统分子动力学模拟。 这一步中，我们将进行整个系统的分子动力学模拟，而不对某些特定原子位置进行限制。因为知识一个小例子，我们将只进行250ps的MD计算。配置文件md2.in如下： -------------------------------------------------------------------------------- oxytocin: 250ps MD &cntrl imin = 0, irest = 1, ntx = 7, ntb = 2, pres0 = 1.0, ntp = 1, taup = 2.0, cut = 10, ntr = 0, ntc = 2, ntf = 2, tempi = 300.0, temp0 = 300.0, ntt = 3, gamma_ln = 1.0, nstlim = 125000, dt = 0.002, ntpr = 100, ntwx = 500, ntwr = 1000 / -------------------------------------------------------------------------------- 上面一些参数解释如下: ntb=2：表示分子动力学过程的压力常数。 ntp＝1：表示系统动力学过程各向同性。 taup = 2.0：压力缓解时间，单位为ps。 pres＝1：引用1个单位的压强。 使用以下命令进行MD： sander –O –i md2.in –o md2.out –p oxy.top –c oxy_md1.rst – r oxy_md2.rst –x oxy_md2.mdcrd –ref oxy_md1.rst –inf md2.info & 模拟时间较长，在P4单机2.0GHz的上需要7.5个小时。 到 此，模拟全部完成，接下来要对得到的数据进行分析。主要数据文件.out文件，包含系统能量、温度，压力等等；.mdcrd文件，是分子动力学轨迹文件， 可以求系统蛋白的RMSD，回转半径等等。数据分析根据不同的研究目的不同而不同，我们将在后文中进行一些简单的分析。 使用pdbviewer软件 load structure的时候，假设侧链缺失会 有提示错误信息。 直接使用xleap产生该结构的top crd后，再使用ambpdb重新产生structure即可加上侧链

本文引用地址： http://www.sciencenet.cn/m/user_content.aspx?id=289032

How to Create a First Shell Script (zz)

How to Create a First Shell Script

http://www.linfo.org/create_shell_1.html

Shell scripts are short programs that are written in a shell programming language and interpreted by a shell process. They are extremely useful for automating tasks on Linux and other Unix-like operating systems. A shell is a program that provides the traditional, text-only user interface for Unix-like operating systems. Its primary function is to read commands (i.e., instructions) that are typed into a console (i.e., an all-text display mode) or terminal window (i.e., all-text mode window) and then execute (i.e., run) them. The default shell on Linux is the very commonly used and highly versatile bash. A programming language is a precise, artificial language that is used to write computer programs, which are sets of instructions that can be automatically translated (i.e., interpreted or compiled) into a form (i.e., machine language) that is directly understandable by a computer's central processing unit (CPU). A feature of bash and other shells used on Unix-like operating systems is that each contains a built-in programming language, referred to as a shell programming language or shell scripting language, which is used to create shell scripts. Among the advantages of using shell scripts are that they can be very easy to create and that a large number are already available in books and on the Internet for use with or without modification for a wide variety of tasks. Shell scripts are also employed extensively in the default installations of Unix-like operating systems. A First Script The following example, although extremely simple, provides a useful introduction to creating and using shell scripts. The script clears the monitor screen of all previous lines and then writes the text Good morning, world. on it. All that is necessary to create this script is to open a text editor (but not a word processor), such as gedit or vi, and type the following three lines exactly as shown on a new, blank page: #!/bin/bash clear echo "Good morning, world." Alternatively, the above code could be copied from this page and pasted to a blank page opened by the text editor page using the standard keyboard or mouse copy and paste functions. After saving this plain text file, with a file name such as morning (or anything else desired), the script is complete and almost ready to run. Scripts are typically run by typing a dot, a forward slash and the file name (with no spaces in between) and then pressing the ENTER key. Thus, for example, if the above script were saved with the name morning, an attempt could be made to execute it by issuing the following command: ./morning However, the script probably will not run, in which case an error message will appear on the screen such as bash: ./morning: Permission denied. This is because the permissions for the file first have to be set to executable. (By default, the permissions for new files are set to read and write only.) The problem can easily be solved by using the chmod command with its 755 option (which will allow the file creator to read, write and execute the file) while in the same directory as that in which the file is located as follows: chmod 755 morning Now the script is ready to run by typing the following, again while in the same directory, and then pressing the ENTER key: ./morning How It Works The first of the three lines tells the operating system what shell to use to interpret the script and the location (i.e., absolute pathname) of the shell. The shell is bash, which is located in the /bin directory (as are all shells); thus the line contains /bin/bash. This instruction is always preceded by a pound sign and an exclamation mark in order to inform the operating system that it is providing the name and location of the shell (or other scripting language). The second line tells the shell to issue the clear command. This is a very simple command that removes all previous commands and output from the console or terminal window in which the command was issued. The third line tells the shell to write the phrase Good morning, world. on the screen. It uses the echo command, which instructs the shell to repeat whatever follows it. (The quotation marks are not necessary in this case; however, it is good programming practice to use them, and they can make a big difference in more advanced scripts.) In slightly more technical terms, Good morning, world. is an argument (i.e., input data) that is passed to the echo command. As is the case with other commands used in shell scripts, clear and echo can also be used independently of scripts. Thus, for example, typing clear on the screen and pressing the ENTER key would remove all previous commands and output and just leave a command prompt for entering the next command. It Doesn't Work! If the phrase Good morning, world. does not appear at the top of the screen, there are several possible reasons: (1) an error was made in copying the code (such as omitting the word echo), (2) the name used in the command was not exactly the same as that of the file (e.g., there is an extra space or a minor difference in spelling or capitalization), (3) the period and/or forward slash were omitted (or reversed) in the command, (4) a space was inserted after the period or slash, (5) the file is not a plain text file (typically because a word processor was used to create it instead of a text editor), (6) the command was not issued in the same directory as that in which the file is located and (7) the permissions were not changed to execute for the owner (i.e., creator) of the file. It is important to avoid practicing writing and executing scripts as the root (i.e., administrative) user. An improperly written script could damage the operating system, and, in a worst case scenario, it could result in the loss of valuable data and make it necessary to reinstall the entire operating system. For this and other reasons, if an ordinary user account does not yet exist on the computer, one should immediately be created (which can be easily accomplished with a command such as adduser). Experiments There are a number of simple, and instructive, experiments that a curious user could do with the above example before moving on to more complex examples. They consist of revising the code as suggested below, saving the revisions (using either the same file name or a different file name), and then executing them as explained above. (1) One is to try changing some of the wording (for example, changing the third line to echo "Good evening, folks."). (2) Another is to add one or more additional lines to be written to the screen, each beginning with the word echo followed by at least one horizontal space. (3) A third is to leave a blank line between two echo lines. (It will be seen that this will not affect the result; however, a blank line can be created by just typing echo on it and nothing else.) (4) A fourth is to insert some blank horizontal spaces. (Notice that the result will be different depending on whether the blank spaces are inserted before or after the first quotation marks. This tells something about the role of quotation marks in shell scripts.) (5) A fifth is to execute the file from a different directory from that in which it is located. This requires adding the path of the executable script to the beginning of the command name when it is issued (e.g., ./test/morning if the file has been moved to a subdirectory named test). (6) Another experiment would be to add some other command to the script file, such as ps (which shows the processes currently on the system), pwd (which shows the current directory), uname (which provides basic information about a system's software and hardware) or df (which shows disk space usage). (Notice that these and other commands can be used in the script with any appropriate options and/or arguments.)

Linux command: Cat

cat [options] [files]

Read (concatenate) one or more files and print them on standard output. Read standard input if no files are specified or if - is specified as one of the files; input ends with EOF.

You can use the > operator to combine several files into a new file, or >> to append files to an existing file. When appending to an existing file, use Ctrl-D, the end-of-file symbol, to end the session.

Options

-A, --show-all

Same as -vET.

-b, --number-nonblank

Number all nonblank output lines, starting with 1.

-e

Same as -vE.

-E, --show-ends

Print $ at the end of each line.

-n, --number

Number all output lines, starting with 1.

-s, --squeeze-blank

Squeeze down multiple blank lines to one blank line.

-t

Same as -vT.

-T, --show-tabs

Print TAB characters as ^I.

-u

Ignored; retained for Unix compatibility.

-v, --show-nonprinting

Display control and nonprinting characters, with the exception of LINEFEED and TAB.

Create file at terminal. To exit, enter STOP.

Examples

cat ch1 Display a file

cat ch1 ch2 ch3 > all Combine files

cat note5 >> notes Append to a file

cat > temp1 Create file at terminal. To exit, enter EOF (Ctrl-D).

cat > temp2 << STOP Create file at terminal. To exit, enter STOP.