richor的个人博客分享 http://blog.sciencenet.cn/u/richor

博文

如何修改gromacs包中的力场

已有 3159 次阅读 2018-7-7 15:31 |个人分类:分子模拟|系统分类:科研笔记|关键词:学者| gromacs, 力场, 修改

【静电修改】

直接去改top文件即可。

【vdw修改】

由于vdw是nonbond相互作用,不能通过直接修改top文件来达到。

需要在include"**/forcefield.itp"和[moleculetypes]之间加入重新定义的原子类型。


比如修改某原子类型Na离子的vdw相互作用(离子通道蛋白)。若不是所有的该原子类型全部都改,那么需要单独定义一个新的原子类型,进而需要重新定义与该原子相关的ffnonbonded.itpffbonded.itp中的相互作用参数。

一个例子脚本如下:

add_bond.sh

#need file ffbonded.itp kcsA_vdw.top, obtain add.top
lbond=6337
lpair=12298
langle=27727
ldihedral1=38499
ldihedral2=54241
lmap=55111

echo "processing the kcsA_vdw.top file"
awk -v lbond=$lbond -v lpair=$lpair '{if(NR>lbond+1 && NR< lpair-2) print}' kcsA_vdw.top > bond.top
awk -v langle=$langle -v ldihedral1=$ldihedral1 '{if(NR>langle+1 && NR< ldihedral1-2) print}' kcsA_vdw.top >angle.top
awk -v ldihedral1=$ldihedral1 -v ldihedral2=$ldihedral2 '{if(NR>ldihedral1+1 && NR <ldihedral2-2) print }' kcsA_vdw.top > dihedral1.top
awk -v ldihedral2=$ldihedral2 -v lmap=$lmap '{if(NR>ldihedral2+1 && NR <lmap-2) print}' kcsA_vdw.top >dihedral2.top

echo "extract the change lines" 
for cx in 899 900 901 902 906 907 908 909 2368 2369 2370 2371 2375 2376 2377 2378 3837 3838 3839 3840 3844 3845 3846 3847 5306 5307 5308 5309 5313 5314 5315 5316
do
grep $cx bond.top >> bond_change.top
done
for cx in 899 900 901 902 906 907 908 909 2368 2369 2370 2371 2375 2376 2377 2378 3837 3838 3839 3840 3844 3845 3846 3847 5306 5307 5308 5309 5313 5314 5315 5316
do
grep $cx angle.top >> angle_change.top
done
for cx in 899 900 901 902 906 907 908 909 2368 2369 2370 2371 2375 2376 2377 2378 3837 3838 3839 3840 3844 3845 3846 3847 5306 5307 5308 5309 5313 5314 5315 5316
do
grep $cx dihedral1.top >> dihedral1_change.top
done
for cx in 899 900 901 902 906 907 908 909 2368 2369 2370 2371 2375 2376 2377 2378 3837 3838 3839 3840 3844 3845 3846 3847 5306 5307 5308 5309 5313 5314 5315 5316
do
grep $cx dihedral2.top >> dihedral2_change.top
done
mv bond_change.top bond.top
mv angle_change.top angle.top
mv dihedral1_change.top dihedral2.top
mv dihedral2_change.top dihedral2.top

echo "obtain the index and name of atom, then substitue the numbers with name"
lstart=71
lend=6336
awk -v lstart=$lstart -v lend=$lend '{if(NR>=lstart && NR<=lend) {if($1!=";")  print $1,$2}}' kcsA_vdw.top > refer.top
awk 'BEGIN{printf("sed \47")} {printf("s/ %s / %s /g;",$1,$2)}END{printf("\47 bond.top >temp\nmv temp bond.top")}' refer.top > refer.sh
echo "do the substitution of bond,may need a little time"
bash refer.sh
awk 'BEGIN{printf("sed \47")} {printf("s/ %s / %s /g;",$1,$2)}END{printf("\47 angle.top >temp\nmv temp angle.top")}' refer.top > refer.sh
echo "do the substitution of angle,may need a little time"
bash refer.sh
awk 'BEGIN{printf("sed \47")} {printf("s/ %s / %s /g;",$1,$2)}END{printf("\47 dihedral1.top >temp\nmv temp dihedral1.top")}' refer.top > refer.sh
echo "do the substitution of dihedral1,may need a little time"
bash refer.sh
awk 'BEGIN{printf("sed \47")} {printf("s/ %s / %s /g;",$1,$2)}END{printf("\47 dihedral2.top >temp\nmv temp dihedral2.top")}' refer.top > refer.sh
echo "do the substitution of dihedral2,may need a little time"
bash refer.sh
rm refer.top
rm refer.sh
echo "finish the substituion"

awk '/CX|HX/{printf("%s %s ",$1,$2);gsub(/CX/,"CA");gsub(/HX/,"HP");printf("%s %s\n",$1,$2)}' bond.top >temp
sort -n temp |uniq > bond.top
awk '/CX|HX/{printf("%s %s %s ",$1,$2,$3);gsub(/CX/,"CA");gsub(/HX/,"HP");printf("%s %s %s\n",$1,$2,$3)}' angle.top >temp
sort -n temp |uniq > angle.top
awk '/CX|HX/{printf("%s %s %s %s ",$1,$2,$3,$4);gsub(/CX/,"CA");gsub(/HX/,"HP");printf("%s %s %s %s\n",$1,$2,$3,$4)}' dihedral1.top >temp
sort -n temp |uniq > dihedral1.top
awk '/CX|HX/{printf("%s %s %s %s ",$1,$2,$3,$4);gsub(/CX/,"CA");gsub(/HX/,"HP");printf("%s %s %s %s\n",$1,$2,$3,$4)}' dihedral2.top >temp
sort -n temp |uniq > dihedral2.top

echo "extracting the lines from ffbonded.itp"
lbondt=1
lconstr=187
langlet=227
ldihed1=701
ldihed2=1303
awk -v lbondt=$lbondt -v lconstr=$lconstr '{if(NR>lbondt+1 && NR<lconstr) print $1,$2,$3,$4,$5}' ffbonded.itp > bond.itp
awk -v lbondt=$lbondt -v lconstr=$lconstr '{if(NR>lbondt+1 && NR<lconstr) print $2,$1,$3,$4,$5}' ffbonded.itp >> bond.itp
sort -n bond.itp | uniq > temp
mv temp bond.itp
awk -v langlet=$langlet -v ldihed1=$ldihed1 '{if(NR>langlet && NR<ldihed1) print $1,$2,$3,$4,$5,$6,$7,$8}' ffbonded.itp > angle.itp
awk -v langlet=$langlet -v ldihed1=$ldihed1 '{if(NR>langlet && NR<ldihed1) print $3,$2,$1,$4,$5,$6,$7,$8}' ffbonded.itp >> angle.itp
sort -n angle.itp | uniq > temp
mv temp angle.itp
awk -v ldihed1=$ldihed1 -v ldihed2=$ldihed2 '{if(NR>ldihed1 && NR<ldihed2) print $1,$2,$3,$4,$5,$6,$7,$8}' ffbonded.itp > dihedral1.itp
awk -v ldihed1=$ldihed1 -v ldihed2=$ldihed2 '{if(NR>ldihed1 && NR<ldihed2) print $4,$3,$2,$1,$5,$6,$7,$8}' ffbonded.itp >> dihedral1.itp
sort -n dihedral1.itp | uniq > temp
mv temp dihedral1.itp
awk -v ldihed2=$ldihed2 '{if(NR>ldihed2+1) print $1,$2,$3,$4,$5,$6,$7}' ffbonded.itp > dihedral2.itp 
awk -v ldihed2=$ldihed2 '{if(NR>ldihed2+1) print $4,$3,$2,$1,$5,$6,$7}' ffbonded.itp >> dihedral2.itp 
sort -n dihedral2.itp | uniq > temp
mv temp dihedral2.itp
echo "processing bondtypes"
echo "[ bondtypes ]" > add.top
wc -l bond.top > lt.txt
lt=`awk '{print $1}' lt.txt`

i=0
until [ $i -eq $lt ]; do
pair=`awk -v i=$i '{if(NR==i+1) printf("%s %s",$3,$4)}' bond.top`
grep "$pair" bond.itp > temp
awk -v i=$i '{if(NR==i+1) print $1,$2}' bond.top > temp1
paste temp temp1 > to_add
awk '{print $6,$7,$3,$4,$5}' to_add >> add.top 
i=$(($i+1))
done
rm bond.top
rm bond.itp

echo "processing angletypes"
echo "[ angletypes ]" >> add.top
wc -l angle.top > lt.txt
lt=`awk '{print $1}' lt.txt`
i=0
until [ $i -eq $lt ]; do
pair=`awk -v i=$i '{if(NR==i+1) printf("%s %s %s",$4,$5,$6)}' angle.top`
grep "$pair" angle.itp > temp
pair1=`cat temp`
if [ -n "$pair1" ]; then
awk -v i=$i '{if(NR==i+1) print $1,$2,$3}' angle.top > temp1
paste temp temp1 > to_add
awk '{print $9,$10,$11,$4,$5,$6,$7,$8}' to_add >> add.top
fi
i=$(($i+1))
done
rm angle.top
rm angle.itp
echo "processing dihedraltypes"
echo "[ dihedraltypes ]" >> add.top
wc -l dihedral1.top > lt.txt
lt=`awk '{print $1}' lt.txt`
i=0
until [ $i -eq $lt ]; do
pair=`awk -v i=$i '{if(NR==i+1) printf("%s %s %s %s",$5,$6,$7,$8)}' dihedral1.top`
grep "$pair" dihedral1.itp > temp
pair1=`cat temp`
if [ -n "$pair1" ]; then
awk -v i=$i '{if(NR==i+1) print $1,$2,$3,$4}' dihedral1.top > temp1
paste temp temp1 > to_add
awk '{print $9,$10,$11,$12,$5,$6,$7,$8}' to_add >> add.top
fi
pair=`awk -v i=$i '{if(NR==i+1) printf("X %s %s X",$6,$7)}' dihedral1.top`
grep "$pair" dihedral1.itp > temp
pair1=`cat temp`
if [ -n "$pair1" ]; then
awk -v i=$i '{if(NR==i+1) printf("X %s %s X\n",$2,$3)}' dihedral1.top > temp1
paste temp temp1 > to_add
awk '{print $9,$10,$11,$12,$5,$6,$7,$8}' to_add >> add.top
fi
i=$(($i+1))
done
rm dihedral1.top
rm dihedral1.itp

echo "[ dihedraltypes ]" >> add.top
wc -l dihedral2.top > lt.txt
lt=`awk '{print $1}' lt.txt`
i=0
until [ $i -eq $lt ]; do
pair=`awk -v i=$i '{if(NR==i+1) printf("%s %s %s %s",$5,$6,$7,$8)}' dihedral2.top`
grep "$pair" dihedral2.itp > temp
pair1=`cat temp`
if [ -n "$pair1" ]; then
awk -v i=$i '{if(NR==i+1) print $1,$2,$3,$4}' dihedral2.top > temp1
paste temp temp1 > to_add
awk '{print $8,$9,$10,$11,$5,$6,$7}' to_add >> add.top
fi
pair=`awk -v i=$i '{if(NR==i+1) printf("X %s %s X",$6,$7)}' dihedral2.top`
grep "$pair" dihedral2.itp > temp
pair1=`cat temp`
if [ -n "$pair1" ]; then
awk -v i=$i '{if(NR==i+1) printf("X %s %s X\n",$2,$3)}' dihedral2.top > temp1
paste temp temp1 > to_add
awk '{print $8,$9,$10,$11,$5,$6,$7}' to_add >> add.top
fi
i=$(($i+1))
done
rm dihedral2.top
rm dihedral2.itp
rm temp1
rm lt.txt

 这里面包含了awk,grep,sed以及bash的if, until, for语句,可以用来参考。

for循环的使用,完全可以替换掉until:

http://gaodr.blog.163.com/blog/static/10461500820093793933887/

建议使用类C的语法:

for ((i=0; i<5; i++))
do
...
done




https://m.sciencenet.cn/blog-485752-1122729.html

上一篇:细菌离子通道的作用
下一篇:关于DNA聚合酶工作机制

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-5-29 01:31

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部