votca-csg 1.3.0-3+b1 / usr / bin / multi_g_rdf

#! /bin/bash
#
# Copyright 2009-2011 The VOTCA Development Team (http://www.votca.org)
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#

g_cmd=${0##*/}
g_cmd=${g_cmd#multi_}
fbase=${g_cmd#g_}
quiet="no"
begin=0
end=""
np=8
index="index.ndx"
out="${fbase}.xvg"
version="1.3 "
single="no"
debug="no"

die() {
  echo -e "$*"
  exit 1
}

qecho() {
  [ "$quiet" = "yes" ] || echo -e "$*"
}

decho() {
  [ "$debug" = "yes" ] && echo -e "$*"
}

show_help () {
  cat << eof
==================================================
========   VOTCA (http://www.votca.org)   ========
==================================================

please submit bugs to bugs@votca.org

${0##*/}, version ${version}



This is a multiplexed version of $g_cmd. (Obsolete will be removed in the next version)

Usage: ${0##*/} [OPTIONS] -- [${g_cmd}_options]

Allowed options:
-N, --NN            Number of tasks
                    Default: $np
-b TIME             Begin time
                    Default: $begin
-e TIME             End time
-n FILE             Name of the index file
                    Default: $index
-o FILE.xvg         Name of the total output file
                    Default: $out
--soutput FILE.xvg  Name of the single output files
                    Default: ${out%.xvg}_NP.xvg (used trunc of name given by -o)
		    (where NP is replaced later by the number of the process)
    --log FILE      Name of logfile
                    Default: ${out%.xvg}_NP.log" (used trunc of name given by -o)
		    (where NP is replaced later by the number of the process)
    --cmd CMD       Change the gromacs command to run
                    Default: ${g_cmd}
    --single        Run only one task at the time
    --debug         Enable debug output
-q, --quiet         Be a little bit quiet
-h, --help          Show this help

Examples:
* ${0##*/} -e 1
* ${0##*/} -e 1 -- -bin 0.05
eof
}

while [ "${1#-}" != "$1" ]; do
 if [ "${1#--}" = "$1" ] && [ -n "${1:2}" ]; then
    #short opt with arguments here: fc
    if [ "${1#-[ebo]}" != "${1}" ]; then
       set -- "${1:0:2}" "${1:2}" "${@:2}"
    else
       set -- "${1:0:2}" "-${1:2}" "${@:2}"
    fi
 fi
 case $1 in
   -[0-9]* | --[0-9]**)
    np=${1#-}
    np=${np#-}
    [ -z "${np//[0-9]}" ] || die "$1 could not be converted to a number"
    shift ;;
   -b)
    begin="$2"
    shift 2;;
   -o)
    out="$2"
    [ -z "${out%%*.xvg}" ] || die "Outfile has to end with xvg"
    shift 2;;
   -n)
    index="$2"
    shift 2;;
   -e)
    end="$2"
    shift 2;;
   --soutput)
    sout="$2"
    [ -z "${sout%%*.xvg}" ] || die "Outfile has to end with xvg"
    shift 2;;
   -q | --quiet)
    quiet="yes"
    shift ;;
   --single)
    single="yes"
    shift ;;
   --cmd)
    g_cmd="$2"
    shift 2;;
   --debug)
    debug="yes"
    shift ;;
   -h | --help)
    show_help
    exit 0;;
   --)
    shift
    break;;
  *)
   die "Unknown option '$1'";;
 esac
done

#defaults
[ -z "$sout" ] && sout="${out%.xvg}_NP.xvg"
[ -z "$logfile" ] && logfile="${out%.xvg}_NP.log"

[ -z "$end" ] && die "-e option is missing!\nHelp with --help"
for i in $g_cmd paste awk sed; do
  [  -n "$(type $i)" ] || die "$i not found"
done

dt="$(awk -v b="$begin" -v e="$end" -v np="$np" 'BEGIN{print (e-b)/np}')"
[ -z "$dt" ] && die "Could not calculate dt"

[ -f "$index" ] || die "Indexfile $index not found"

cat <<EOF
${0##*/} is obsolete. Please use csg_stat instead. We will remove ${0##*/}
in the next version of VOTCA. csg_stat is parallel as well (--nt option)
and uses a smarter neighbor search algrithm (grid search), which speeds
the caluclation up a lot. Additionally it supports calculations for
multiple interaction (incl. bonded) at once. The only thing you have
to do, is to create as xml file like this:

<cg>
  <!-- example for a non-bonded interaction entry -->
  <non-bonded>
    <!-- name of the interaction -->
    <name>A-A</name>
    <min>0.0</min>
    <max>2.23</max>
    <step>0.01</step>
    <!-- these are atom TYPES not energy groups, * (e.g. A*) is allowed. -->
    <type1>A1</type1>
    <type2>A1</type2>
  </non-bonded>
  <!-- another non-bonded block for another interaction -->
  <non-bonded>
    ..
  </non-bonded>
</cg>

and execute:
csg_stat -top topol.tpr --trj traj.trr --nt XXX --options FILE.xml
EOF
sleep 30

echo "Select a reference group and 1 group"
sed -n 's/.*\[\(.*\)\].*/\1/p' $index | awk '{print "Group",NR-1,"(",$0,")"}'
echo -n "Select a group: "
read -r || die "Wrong input"
input="${REPLY}"
echo -n "Select a group: "
read -r || die "Wrong input"
input="${input}\n${REPLY}"
qecho ""
qecho "Selected: $input"

b="$begin"
for ((i=0;i<$np;i++)); do
  np0="$(printf "%02i" $i)"
  this_log="${logfile/NP/$np0}"
  this_out="${sout/NP/$np0}"
  e="$(awk -v b="$begin" -v dt="$dt" -v i="$i" 'BEGIN{print dt*(i+1)+b}')"
  echo -e "$input" | $g_cmd -b $b -e $e -n $index -o $this_out $* &> $this_log &
  pid="$!"
  qecho "Calculation of $fbase from $b ps to $e ps started (pid $pid, logfile $this_log)"
  b="$e"
  if [ "$single" = "yes" ]; then
    wait $pid || die "The $fbase calculation with pid $failed failed"
  else
    pids="${pids} $pid"
  fi
done

if [ "$single" = "no" ]; then
  qecho "Waiting for single $fbase calculations to finish"
  failed=""
  for i in $pids; do
    wait $i || failed="$failed $i"
  done
  [ -z "$failed" ] || die "The $fbase calculations with pids $failed failed"
fi

qecho "Checking output"
lines=""
for ((i=0;i<$np;i++)); do
  np0="$(printf "%02i" $i)"
  this_out="${sout/NP/$np0}"
  [ -f "$this_out" ] || die "Could not find $this_out"
  this_lines="$(sed -n '$=' "$this_out")" || die "linecount of $this_out failed"
  [ -z "$this_lines" ] && die "linecount of $this_out was empty"
  [ -z "$lines" ] && lines="$this_lines"
  [ $this_lines -eq $lines ] || die "Number of data lines differs for $fbase nr $i - $this_lines vs $lines lines"
done

qecho "Calculating average"
pasteline=
mean="0"
mean2="0"
mco="\$$((2*$np+1))"
for ((i=0;i<$np;i++)); do
  np0="$(printf "%02i" $i)"
  this_out="${sout/NP/$np0}"
  co="\$$((2 * $i + 2 ))"
  pasteline="$pasteline $this_out"
  mean="${mean}+${co}"
  mean2="${mean2}+($co-$mco)**2"
done
awkprog="/[@#]/{next;}{print \$0,($mean)/$np}"
awkprog2="/[@#]/{next;}{print \$1,$mco,sqrt(($mean2)/$(($np*($np-1))))}"
decho "awk code to calc mean: $awkprog"
decho "awk code to calc block error: $awkprog2"
paste $pasteline | awk "$awkprog" | awk "$awkprog2" > $out