GROMACS自动化SHELL操作基础

基本GROMACS的SHELL操作方式

其实这个在GROMACS官网有一些介绍Using Commands in Scripts,但是命令比较旧,可能不太全面明白。我在这里简单介绍一些方法。

一、比如我想计算蛋白组名(Protein),编号(1),那么我们可以这样

    1. echo方式

      1
      echo 1 1 |gmx rms -f md.trr -s md.tpr

      或者

      1
      echo Protein Protein |gmx rms -f md.trr -s md.tpr
    1. <<EOF方式
      1
      2
      3
      4
      gmx rms -f md.trr -s md.tpr <<EOF
      1
      1
      EOF
    1. 文本输入重定向方式
      建立一个rmsd.txt文本文件,内容如下:

      1
      2
      1
      1

      然后输入命令:

      1
      gmx rms -f md.trr -s md.tpr < rmsd.txt

二、需要以Ctrl+D结束的命令,比如计算水分子组名(Water),编号(12)的rdf

    1. echo方式

      1
      echo -e "12\n12" |gmx rdf -f md.trr -s md.tpr

      或者

      1
      echo -e "Water\nWater" |gmx rdf -f md.trr -s md.tpr

      这里的\n相当于你手动敲击回车。

    1. 直接通过工具的选择命令:
      1
      gmx rdf -f md.trr -s md.tpr -ref "group 12" -sel "group 12"
    1. <<EOF方式
      1
      2
      3
      4
      gmx rdf -f md.trr -s md.tpr <<EOF
      12
      12
      EOF
    1. 文本输入重定向方式
      建立一个rdf_water.txt文本文件,内容如下:

      1
      2
      12
      12

      然后输入命令:

      1
      gmx rdf -f md.trr -s md.tpr < rdf_water.txt

SHELL之并行

以博文使用GROMACS计算径向分布函数为例。

这篇博文讲解中包括了一个并行部分,这也是在数据处理中十分常见的。我详细解读了一番,发现还有一定的扩展空间。SHELL调用外部程序,如果本身程序不支持并行,那么我们有几种方式来进行任务并行:

  • 如上面博文所述,使用paralle。但是这个程序需要你自己额外安装,因此不是十分方便。
  • 直接丢任务到后台,即在循环运行命令后面加&即可。这种方式十分的危险,特别是任务耗时长,CPU核心数少的情况。这样会很容易导致后台任务过多,使得计算机瘫痪!
  • 将任务分到每一个CPU核心上,即满载,等待所有核心完成任务再进行下一批任务。这种方法的实现很简单,只需要确认机器有多少核心即可,脚本如下:
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    #!/bin/bash
    startTime=`date +%Y%m%d-%H:%M:%S`
    startTime_s=`date +%s`
    nprocs=$(cat /proc/cpuinfo |awk '/^processor/ {print $3}'|wc -l)

    mkdir rdf
    i=0
    for j in {3..500}
    do
    #echo "Computing RDF between group $j - group 2"
    echo -e "$j\n2" |gmx rdf -f md.xtc -s md.tpr -n index.ndx -o ./rdf/rdf_${j}_2.xvg -quiet > /dev/null 2>&1 &
    let i+=1
    [[ $((i%nprocs)) -eq 0 ]] && wait
    done

    echo -e "\nWaiting all jobs finish..."
    wait
    echo "All jobs have finished!"

    endTime=`date +%Y%m%d-%H:%M:%S`
    endTime_s=`date +%s`
    sumTime=$[ $endTime_s - $startTime_s ]
    echo "$startTime ---> $endTime" "Totl:$sumTime seconds"

上面程序通过判断每次分摊到CPU核上的任务是否全部完成,然后进行下一批。此方式适合每个任务耗时相差不大的时候使用,如果耗时相差很大,就会导致任务堵塞,慢得很!

  • 时刻检查空闲CPU内核,将任务实时传给空闲核心。这个看上去就比较的复杂了,但是效率更高,参考Parallel batch processing in the shell。代码如下:
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    #!/bin/bash
    startTime=`date +%Y%m%d-%H:%M:%S`
    startTime_s=`date +%s`

    function mapp() {
    if [[ -z $MAPP_NR_CPUS ]] ; then
    local MAPP_NR_CPUS=$(cat /proc/cpuinfo |awk '/^processor/ {print $3}'|wc -l)
    fi
    local mapp_pid=$(exec bash -c 'echo $PPID')
    local mapp_funname=$1
    local -a mapp_params
    mapp_params=("$@")
    local mapp_nr_args=${#mapp_params[@]}
    local mapp_current=0
    function mapp_trap() {
    echo "MAPP PROGRESS: $((mapp_current*100/mapp_nr_args))%" 1>&2
    if [[ $mapp_current -lt $mapp_nr_args ]] ; then
    let mapp_current+=1
    (
    $mapp_funname "${mapp_params[$mapp_current]}"
    kill -USR1 $mapp_pid
    ) &
    fi
    }

    trap mapp_trap SIGUSR1
    while [[ $mapp_current -lt $mapp_nr_args ]]; do
    wait
    if [[ $mapp_current -lt $mapp_nr_args && $? -lt 127 ]] ; then
    sleep 1
    local mapp_tmp_count=$mapp_current
    wait
    if [[ $mapp_tmp_count -eq $mapp_current ]] ; then
    echo " MAPP_FORCE" 1>&2
    for i in $(seq 1 ${MAPP_NR_CPUS}) ; do
    mapp_trap
    done
    fi
    fi
    done
    for i in $(seq 1 ${MAPP_NR_CPUS}) ; do
    wait
    done
    trap - SIGUSR1
    unset -f mapp_trap
    }

    function rdf() {
    #echo "Computing RDF between group $j - group 2"
    echo -e "$1\n2" |gmx rdf -f md.xtc -s md.tpr -n index.ndx -o ./rdf/rdf_${1}_2.xvg -quiet > /dev/null 2>&1
    }

    mkdir rdf
    (mapp rdf {3..500})

    endTime=`date +%Y%m%d-%H:%M:%S`
    endTime_s=`date +%s`
    sumTime=$[ $endTime_s - $startTime_s ]
    echo -e "\n$startTime ---> $endTime" "Totl: $sumTime seconds"

此代码适合每个任务耗时不一样的时候使用,比前一个代码更具有效率