localize.sh

#!/bin/bash
'''Script taken from CamCASP 5.5 software, and original author contributions
are as listed in 
http://www-stone.ch.cam.ac.uk/programs/camcasp.html
'''

function help {
cat <<EOF
$0
Uses the Lillestolen-Wheatley and pfit procedures to localize a set of
polarizability files and refine the local polarizabilities.

Usage:
  ${0##/*/} [--name] <name> [options]
or
  ${0##/*/} -i <infile> --pol <pol-file-prefix> -o <output pol-file-prefix> \
      [options]

The first form is sufficient for files prepared using the standard
procedure. In this case the values are set to infile=<name>.ornt,
input pol-file-prefix=<name>_NL4_ and output pol-file-prefix=<name>_L${limit}_ 
respectively, where limit is the specified maximum rank (default 2).

Options:
  --axes <file>
        The local axes are defined in the specified file. Default file
        is <name>.axes; if this file doesn't exist in the current directory
        and there is a <name>.axes file in the directory above, that is
        used. If <name>.axes doesn't exist in either of these places or is
        empty, global axes are used.
        Note that global axes must be used for the localization step; the
        local axes, if any, are applied to the final polarizabilities and
        dispersion coefficients.
  --sites <file>
        The molecular sites are defined in the specified file. Default file
        is name.sites.
  --keep
        Don't delete the intermediate files generated by the script.
        Shouldn't be necessary unless bugs occur.
  --limit <limit>
        Generate local polarizabilities up to rank <limit>. Default 2.
  --wsmlimit <n>
        Limit the rank of refined local polarizabilities to n. Must be
        less than or equal to <limit>. Default same as <limit>.
  --hlimit <n>
        On hydrogen sites, generate refined local polarizabilites up to
        rank n. <hlimit> must not be greater than <wsmlimit>. Default
        same as <wsmlimit>.
  --weight <w>
        Use penalty weighting scheme <w> for refining the polarizabilities.
        See the CamCASP User's Guide for details.
  --isotropic
        When refining, fit only isotropic polarizabilities. Default is
        to fit anisotropic terms as well.
  --norefine
        Skip the refinement step.
  --nodisp
        Skip the calculation of the dispersion coefficients.
  --loc
        Localization method. Either LW (Lillestolen-Wheatley, default)
        or LS (LeSueur-Stone).
  --svd
        Use singular value decomposition, discarding singular values
        smaller than 1D-5*largest.

Steps:

In the following, the file names in parentheses are those constructed if
only <name> is specified.

The script looks for a .pol file in the current directory or in the
OUT subdirectory. If there is more than one, the user is asked to
choose between them. Similarly it looks for a .p2p file containing the
point-to-point polarizabilities.

The polarizability file is split into sections, one for each
frequency. It is assumed that the static polarizabilities are
included, followed by dynamic polarizabilities at 10 frequencies.
The file names are <prefix>nnn.pol, with nnn = 000, 001, ..., 010.

The orient input file <infile> (<name>.ornt) is edited to specify the 
local-axes file and the sequence number, and Orient is run to carry out 
the localization at each frequency.

These single-frequency files are concatenated into a single file
<output pol-file-prefix>_0f10.pol (<name>_L<limit>_0f10.pol).

Now the point-to-point polarizabilities are split into
single-frequency files, and the PFIT program is run for each one, to
refine the local polarizabilities.

Finally these refined single-frequency files are concatenated into a
single file.

EOF
}

[[ $# -eq 0 || "$1" = "--help" ]] && help && exit

. $CAMCASP/bin/functions

bin=$CAMCASP/bin
  

file=""
name=""
axes=""
sites=""
limit=2
hlimit="1"
wsmlimit=""
isotropic=""
infile=""
prefix=""
keep=""
out=""
disp=""
refine=""
weight=4
loc="LW"
svd="OFF"
while [[ $# -gt 0 ]]; do
  case $1 in
    --name )
      name=$2
      shift;;
    --axes )
      axes=$2
      shift;;
    --sites )
      sites=$2
      shift;;
    --limit )
      limit=$2
      shift;;
    --hlimit )
      hlimit=$2
      shift;;
    --wsmlimit )
      wsmlimit=$2
      shift;;
    --weight )
      weight=$2
      shift;;
    --isotropic )
      isotropic="ISOTROPIC";;
    -i )
      infile=$2
      shift;;
    -o )
      out=$2
      shift;;
    --pol )
      prefix=$2
      shift;;
    --keep )
      keep="yes";;
    --nodisp )
      disp="no";;
    --norefine )
      refine="no";;
    --loc )
      loc=$2
      shift;;
    --svd )
      svd="threshold 1d-5";;
    --* )
      die "I don't understand $1";;
    * )
      if [[ -z $name ]]; then
        name=$1
      else
        die "I don't understand $1"
      fi;;
  esac
  shift
done

[[ -z $name ]] && die "File prefix (system name) not specified"
[[ -z $wsmlimit ]] && wsmlimit=$limit

[[ $wsmlimit -gt $limit ]] && die "wsmlimit must not be greater than limit"
[[ $hlimit -gt $wsmlimit ]] && die "hlimit must not be greater than wsmlimit"

rm -f error_file error.log pfit.error

PS3="Enter number for required file: "
#  Look for polarizability file
findfile .pol OUT || findfile .pol || die "No .pol file found"
polfile=$found
echo "Using polarizability file $polfile"
  
#  Look for p2p file
findfile .p2p OUT || findfile .p2p || die "No .p2p file found"
p2pfile=$found
echo "Using p2p file $p2pfile"
  
[[ -z $infile ]] && infile=${name}.ornt
[[ -z $prefix ]] && prefix="${name}_NL4_"
[[ -z $out ]] && out="${name}_L${limit}_"
[[ -z $polfile ]] && die "Polarizability file not specified"
[[ -e $infile ]] || die "Orient input file $infile not found"
#  Axis definition file: if not present, link to file in directory above,
#  if present, otherwise create empty file.
[[ -z $axes ]] && axes=$name.axes
[[ ! -e $axes ]] && [[ -e ../$axes ]] && ln -s ../$axes .
[[ -e $axes ]] || touch $axes

echo -n "Splitting $polfile into individual frequencies ..."
echo "/usr/bin/csplit --prefix=$prefix --suffix-format="%03d.pol" \
   --elide-empty-files --quiet $polfile "/# INDEX/" '{*}'"
/usr/bin/csplit --prefix=$prefix --suffix-format="%03d.pol" \
   --elide-empty-files --quiet $polfile "/# INDEX/" '{*}'
echo " done"
cp ${prefix}000.pol ${prefix}static.pol

echo "Localizing polarizabilities using $loc procedure ..."
rm -f ${out}0f10.pol
for tag in 000 001 002 003 004 005 006 007 008 009 010
do
echo -n "$tag "
  if [[ -e $prefix$tag.pol ]]; then
    #  Localize the polarizabilities, up to rank $limit, for this frequency index
    #  Write them to $out$tag.pol
    $bin/replace AXES=$axes PREFIX=$prefix INDEX=$tag LIMIT=$limit \
        LOC=$loc < $infile | orient > $out$tag.out
    [[ -e orient.error ]] && die "Error in localization"
    [[ $tag = "000" ]] && cp -p $out$tag.pol ${out}static.pol
    #  Concatenate the localized files
    cat $out$tag.pol >> ${out}0f10.pol || die "error"
    [[ -z $keep ]] && rm $out$tag.out $out$tag.pol $prefix$tag.pol
  else
    if [[ $tag = "000" ]]; then
      echo "Warning -- no static polarizabilities found"
    else
      die "Can't find polarizability file $prefix$tag.pol"
    fi
  fi
done
echo " ... done"

if [[ -z $refine ]]; then

#  Site definition file: if not present here, link to file in directory
#  above if present there.
[[ -z $sites ]] && sites=$name.sites
[[ -e $sites ]] || ln -s ../$sites . || die "Site definition file $sites not found"
#  Polarizability model definition: if not present here but present in
#  directory above, link to that.
pdef=${name}.pdef
[[ ! -e $pdef ]] && [[ -e ../$pdef ]] && ln -s ../$pdef .

#  Have the p2p files already been unpacked?
ok=""
for tag in 000 001 002 003 004 005 006 007 008 009 010
do
  [[ -e ${name}_${tag}.p2p ]] || ok="no"
done
if [[ -n $ok ]]; then
echo "Splitting the point-to-point polarizability file ..."
#  Split the point-to-point polarizability file
rm -f ${name}_???.p2p
process <<EOF
read P2P pols for $name
  Frequencies Static + 10
  P2P-Pols $p2pfile SPLIT
End
EOF
echo "done"
fi

echo "Refining the local polarizabilities"
rm -f ${name}_ref_wt${weight}_L${wsmlimit}_0f10.pol
for tag in 000 001 002 003 004 005 006 007 008 009 010
do
    echo -n "$tag "
    n=${tag##0}; n=${n##0}  #  Strip leading zeroes
    $bin/replace CAMCASP=$CAMCASP INDEX=$n LIMIT=$limit \
    HLIMIT=$hlimit WSMLIMIT=$wsmlimit ISOTROPIC=$isotropic WEIGHT=$weight \
       < $name.prss > $name.$tag.prss
       process < $name.$tag.prss | \
       $bin/replace AXES=$axes SITES=$sites PDEF=$pdef SVD="$svd" > ${name}_pfit_$tag.data
    [[ -e error_file || -e error.log ]] && die "Error in process"
    # Following line now redundant.
    # [[ -e ${name}.pdef ]] && $bin/replace_poldefs.pl ${name} $tag
    pfit < ${name}_pfit_$tag.data > ${name}_ref_wt${weight}_L${wsmlimit}_$tag.out
    [[ -e pfit.error ]] && die "Error in pfit -- see file pfit.error"
    echo "# INDEX $tag" >> ${name}_ref_wt${weight}_L${wsmlimit}_0f10.pol
    cat ${name}_ref_wt${weight}_L${wsmlimit}_f$n.pol >> ${name}_ref_wt${weight}_L${wsmlimit}_0f10.pol
    echo " " >> ${name}_ref_wt${weight}_L${wsmlimit}_0f10.pol
    if [[ $n = "0" ]]; then
      echo "# Static polarizabilities" > ${name}_ref_wt${weight}_L${wsmlimit}_static.pol
      cat ${name}_ref_wt${weight}_L${wsmlimit}_f0.pol >> ${name}_ref_wt${weight}_L${wsmlimit}_static.pol
    fi
    [[ -z $keep ]] && rm ${name}_$tag.p2p ${name}_pfit_$tag.data* \
      ${name}_ref_wt${weight}_L${wsmlimit}_f$n.pol ${name}.$tag.prss
done
echo " ... done"
else
  echo "Refinement omitted"
fi

if [[ -z $disp && -z $refine ]]; then
    echo "Calculating the dispersion coefficients ..."
    [[ -e casimir.error ]] && rm casimir.error
    casimir_in=${name}_ref_wt${weight}_L${wsmlimit}.casimir
    casimir_out=${name}_ref_wt${weight}_L${wsmlimit}_casimir.out
    case $wsmlimit in
        "1" ) maxN="6"  ;;
        "2" ) maxN="10" ;;
        "3" ) maxN="12" ;;
          * ) maxN="n"  ;;
    esac
    if [[ $isotropic == "ISOTROPIC" ]]; then
      potfile="${name}_wt${weight}_L${wsmlimit}iso_C${maxN}iso.pot"
    else
      potfile="${name}_wt${weight}_L${wsmlimit}_C${maxN}.pot"
    fi
    $bin/replace WSMLIMIT=$wsmlimit HLIMIT=$hlimit WEIGHT=$weight < ${name}_casimir.prss | \
    process > ${casimir_in}
    casimir < ${casimir_in} > ${casimir_out}
    [[ -e casimir.error ]] && die "Error in casimir -- see file casimir.error"
    #Now create the potential file
    perl -e "while (<>) {last if /Dispersion/} while (<>) {print}" \
      < ${casimir_out} > ${potfile}
    [[ -s ${potfile} ]] || die "Dispersion coefficient calculation failed"
    # csplit --prefix="Cn_disp" --suffix-format="%1d.coeff" ${casimir_out} \
    #    "%Dispersion%+1" || die "Dispersion coefficient calculation failed"
    # mv Cn_disp0.coeff ${potfile}
    #
    echo "...done."
    echo "Dispersion coefficients are in ${casimir_out}"
    echo "The dispersion potential, in Orient form, is in ${potfile}"
fi

exit