计算化学公社

 找回密码 Forget password
 注册 Register
Views: 11519|回复 Reply: 2
打印 Print 上一主题 Last thread 下一主题 Next thread

[GROMACS] 求助:重复Gromacs官网上的伞形抽样教程得到的summary_distances.dat文件出现问题

[复制链接 Copy URL]

16

帖子

0

威望

1920

eV
积分
1936

Level 5 (御坂)

跳转到指定楼层 Go to specific reply
楼主
大家好,我得到的summary_distances.dat文件中的距离不对,而且还不连续,到最后就没有了,谢谢!

summary_distances.dat

10.16 KB, 下载次数 Times of downloads: 17

313

帖子

2

威望

3900

eV
积分
4253

Level 6 (一方通行)

2#
发表于 Post on 2021-3-29 19:41:34 | 只看该作者 Only view this author
本帖最后由 lyj714 于 2021-3-29 19:45 编辑

如果你用的是最原始的那个perl脚本,那是因为以前的这个脚本有些问题。如果是用的现在原英文教程中的那个bash脚本就没啥问题,不过速度慢。本人提供一个我自用的脚本,可以一步得到US的批处理脚本
  1. #!/bin/bash
  2. # A awk script for creating batch script of Umbrella Sample, GMX Version should >= 5.1

  3. #vvvvvvvvvvvvvvvvvvvvvvvvvvvCHANGE SECTIONvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv#
  4. #From which step to process, support 1 and 2, 1 represents need tpr, xtc/trr and ndx file
  5. #not need summary_distances.dat file; 2 represents only need summary_distances.dat file
  6. # summary_distances.dat includs:
  7. # two columns data
  8. #      1st column: frame number, from gmx trjconv -sep
  9. #      2nd column: distance [nm] or angle [degree], from gmx distance or gmx gangle
  10. step=1                                  # starting step

  11. # For example, tpr and xtc/trr file name, pull group name
  12. # for calculating center of mass distance/dihedral should be changed by yourself
  13. tpr="pull.tpr"                          # pull tpr file
  14. xtc="pull.xtc"                          # pull xtc/trr file
  15. groupA="Chain_A"                        # first group name for distance pull or dihedral pull(noly)
  16. groupB="Chain_B"                        # second group name for distance pull

  17. # gmx name, interval, ndx, pull type
  18. # if step=2, only need input file [summary_distances.dat], also can see gromacs tutorials
  19. gmx="gmx"                               # Maybe gmx, gmx_mpi etc
  20. inputfile="summary_distances.dat"       # temp file, do not change
  21. ndx="index.ndx"                         # ndx file, including groupA, groupB
  22. pull_type="distance"                    # distance pull, maybe "dihedral"
  23. interval=0.2                            # Unit nm for distance, degree for dihedral, you can change it
  24. outdist="outdistance.dat"               # distance/dihedral output file
  25. outbatch="do_batch.bsh"                 # batch output file
  26. #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#

  27. if [[ $step -le 1 ]];then
  28. #============================DO NOT CHANGE==============================#
  29. # generate conf*.gro file in order to US simulation
  30. echo -e "\n>>Generating a series of conf*.gro file..."
  31. echo 0 | $gmx trjconv -s $tpr -f $xtc -o conf.gro -sep -quiet > /dev/null 2>&1

  32. # generate distance/dihedral of each frame and change the first column into frame index(from 0 start)
  33. if [[ $pull_type"x" == "distance""x" ]]; then
  34.     echo -e "\n>>Calculating center of mass distance between two groups..."
  35.     $gmx distance -s $tpr -f $xtc -n $ndx -select "com of group "$groupA" plus com of group "$groupB"" -oall temp.xvg -xvg none -quiet > /dev/null 2>&1
  36. elif [[ $pull_type"x" == "dihedral""x" ]];then
  37.     echo -e "\n>>Calculating dihedral ..."
  38.     echo $groupA |$gmx gangle -s $tpr -f $xtc -n $ndx  -g1 dihedral -oall temp.xvg -xvg none -quiet > /dev/null 2>&1
  39. fi
  40. awk '{printf("%4d %8f\n", NR-1, $2)}' temp.xvg > $inputfile
  41. fi

  42. if [[ $step -le 2 ]];then
  43. # read input file and generate distance/dihedral file and batch file for simulation
  44. echo -e "\n>>Creating batch script according to given interval..."
  45. awk '
  46.     BEGIN {
  47.         inputfile = "'$inputfile'"; interval = "'$interval'"+0;
  48.         outdist   = "'$outdist'";   outbatch = "'$outbatch'";
  49.         gmx = "'$gmx'"; ndx = "'$ndx'"; pull_type = "'$pull_type'"
  50.         print "\n>>Input file:     ", inputfile
  51.         print ">>Input interval: ", interval, "(nm)"

  52.         #check file
  53.         Check_File(inputfile)

  54.         #read file and obtain data
  55.         ll = sampledist(inputfile, interval, sample_indexs, distances)

  56.         #output data file, including distance and batch script
  57.         if(pull_type == "distance")
  58.             printf("%10s %10s %10s\n", "frame", "dist", "d_dist")   > outdist
  59.         else if(pull_type == "dihedral")
  60.             printf("%10s %10s %10s\n", "frame", "angle", "d_angle") > outdist
  61.         printf("#!/bin/bash\n\n")                                   > outbatch
  62.         for(i = 0; i < ll; i++) {
  63.             frame = distances[sample_indexs[i], 0]
  64.             dist  = distances[sample_indexs[i], 1]
  65.             if(i == 0) {delta_dist = "NA"}
  66.             else {
  67.                 prev_dist  = distances[sample_indexs[i-1], 1]
  68.                 delta_dist = dist - prev_dist
  69.             }

  70.             out_dist(frame, dist, delta_dist, outdist)
  71.             out_batch(frame, outbatch)
  72.         }

  73.         Thanks()
  74.     }

  75.     #==========================SOME FUNCTION============================#
  76.     #output batch file, you can change it if you understand it means
  77.     function out_batch(fram, outfile) {
  78.         printf("#Do the %s frame simulation\n", fram)                   > outfile
  79.         printf("%s grompp -f npt_umbrella.mdp -c conf%s.gro -r conf%s.gro -p topol.top -n %s -o npt%s.tpr\n",
  80.                 gmx, fram, fram, ndx, fram)                             > outfile
  81.         printf("%s mdrun -deffnm npt%s\n", gmx, fram)                   > outfile
  82.         printf("%s grompp -f md_umbrella.mdp -c npt%s.gro -r npt%s.gro -t npt%s.cpt -p topol.top -n %s -o umbrella%s.tpr\n",
  83.                 gmx, fram, fram, fram, ndx, fram)                       > outfile
  84.         printf("%s mdrun -deffnm umbrella%s\n\n", gmx, fram)            > outfile
  85.     }

  86.     #output distance file
  87.     function out_dist(fram, dist, delta_dist, outdist) {
  88.         printf("%10d %10.3f %10s\n", fram, dist, delta_dist)            > outdist
  89.     }

  90.     #main function for finding appropriate windows
  91.     function sampledist(table, interval, sample_indexs, distances, _ARGVEND_, i, n) {
  92.         n = 0;
  93.         while(getline<table>0) {
  94.             distances[n, 0] = $1
  95.             distances[n, 1] = $2+0
  96.             n++
  97.         }
  98.         close(table)

  99.         current_index = 0
  100.         sample_indexs[0] = 0
  101.         iindex = 0
  102.         while(current_index < n) {
  103.             target_distance = distances[current_index, 1] + interval
  104.             tt = 0
  105.             for(i = current_index; i < n; i++) {
  106.                 onward[tt++]  = abs(target_distance - distances[i, 1])
  107.             }
  108.             mim_value = min(onward)
  109.             for(j in onward) {
  110.                 if(onward[j] == mim_value) {
  111.                     next_index = j + current_index
  112.                     break
  113.                 }
  114.             }
  115.             if(current_index == next_index) {break}
  116.             else {
  117.                 sample_indexs[++iindex] = next_index
  118.                 current_index = next_index
  119.             }
  120.             delete onward
  121.         }
  122.         return length(sample_indexs)
  123.     }

  124.     #abs function
  125.     function abs(a) { return a>0? a : -a }

  126.     #find minimum
  127.     function min(x, _ARGVEND_, value, i) {
  128.         value = 99999
  129.         for(i in x) {
  130.             value = value < x[i] ? value : x[i]
  131.         }
  132.         return value
  133.     }

  134.     #output information
  135.     function Thanks() {
  136.         print "\n=============RESULT=============="
  137.         print outdist " has been created!"
  138.         print outbatch " has been created!"
  139.         print "================================="
  140.     }

  141.     # check file
  142.     function Check_File(file)
  143.     {
  144.         if (getline < file < 0)
  145.         {
  146.             printf("\nError! file [%s] not found!", file)
  147.             exit 1
  148.         }
  149.         else
  150.             close(file)
  151.     }
  152. '
  153. fi

复制代码

16

帖子

0

威望

1920

eV
积分
1936

Level 5 (御坂)

3#
 楼主 Author| 发表于 Post on 2021-3-30 00:21:22 | 只看该作者 Only view this author
lyj714 发表于 2021-3-29 19:41
如果你用的是最原始的那个perl脚本,那是因为以前的这个脚本有些问题。如果是用的现在原英文教程中的那个ba ...

收到,非常感谢!

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2026-2-22 17:36 , Processed in 0.166643 second(s), 23 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list