#!/bin/bash
echo "================= cook configurator (part of MACSIMUS) ====================="
echo "|  This script must be called from the MACSIMUS/cook/ directory!            |"
echo "|  The cook version will be configured and compiled in a SUBDIRectory.     |"
echo "|  Synopsis:                                                               |"
echo "|    ./configure.sh [SUBDIR [EXENAME]]                                     |"
echo "|  where EXENAME must not contain spaces and other control characters      |"
echo "|  You will be asked for compile-time options:                             |"
echo "|    Answer by a letter or a number; just <Enter> means \"no\" or \"default\"  |"
echo "============================================================================"

METAMAKE=metamake.mmk

if [ -e ../makemake.sh ] ; then
  source ../makemake.sh
  echo "reading ../makemake.sh"
else
  echo "../makemake.sh not found - MACSIMUS has not been properly configured!"
  echo "Enter compiler name:"
  read
  CC=$REPLY
fi

if [ "$CC" = "" ]; then
  CC=cc
  echo "NOTE: $CC will be passed to makemake"
fi

echo NOTE: compiler CC=$CC

EXENAME=""
if [ $# -gt 0 ] ; then
  SUBDIR="$1"
  echo "NOTE: \"$SUBDIR/\" specified as subdirectory name"
  if [ $# -gt 1 ] ; then
    EXENAME="$2"
  echo "NOTE: \"$EXENAME\" specified as the executable name"
  fi
else
  echo -n "Enter subdirectory name (will be created if necessary): "
  read
  SUBDIR="$REPLY"
  if [ "$SUBDIR" = "" ] ; then
    echo "no directory specified - stop"
    exit
  fi
fi

if [ ! -e "$SUBDIR" ] ; then
  mkdir -p "$SUBDIR"
  echo "NOTE: \"$PWD/$SUBDIR\" created"
else
  if [ -d "$SUBDIR" ] ; then
    echo "NOTE: directory \"$PWD/$SUBDIR/\" already exists"
  else
    echo "ERROR: \"$SUBDIR\" exists and is not a directory"
    exit
  fi
fi

COOKDIR="$PWD"
cd "$SUBDIR"

NAME="cook"
MODULE=""
PROJECT="COOK"

if [ -e simopt.h ] ; then
  mv simopt.h simopt.h~
  echo "NOTE: existing (previous) simopt.h was renamed to simopt.h~"
fi

echo "/* This file was generated by configure.sh */" > simopt.h
echo "/* See cook/generic/simopt.h for more details */" >> simopt.h
echo "#define COOK" >> simopt.h

echo
echo "===== Select BOUNDARY CONDITIONS:"
echo "> p = Periodic in 3D, spherical cut-off of Lennard-Jones interactions"
echo "      Allows later modifications (wall, slab, periodic sum) (DEFAULT)"
echo "  f = Free (vacuum), true electrostatic 1/r (no cutoff)"
echo "      NB: not compatible with Gaussian charges"
echo "  s = Slit, as above but periodic in x,y only; NOT TESTED RECENTLY"
echo "  n = Periodic in 3D, all pair interactions use nearest-image distance"
echo "      (special, cutoff corections not implemented; NOT TESTED RECENTLY)"
read
BC=$REPLY
if [ "$BC" = "" ] ; then
  BC=p
fi

export COULOMB=COULOMB_undefined

case "$BC" in
n )
  NAME=$NAME"ni"
  PROJECT=$PROJECT"-ni"
  MODULE=$MODULE" ewald"
  COULOMB=0 # to ask correctly for widom
  echo "#define NIBC" >> simopt.h
  ;;
f )
  echo
  echo "===== Select charge-charge/external interaction:"
  echo "> c = Coulomb interaction (DEFAULT)"
  echo "  s = Stars (gravity: same charges attract each-other, #define STARS)"
  echo "  g = Gaussian instead of force field (#define GAUSSIAN)"
  echo "      NB: many measurements are missing/out of order"
  read
  FREE=$REPLY
  if [ "$FREE" = "" ] ; then
    FREE="c"
  fi
  echo "#define FREEBC" >> simopt.h
  case "$FREE" in
   c )
      NAME=$NAME"free"
      PROJECT=$PROJECT"-free"
     ;;
   s )
      NAME=$NAME"stars"
      PROJECT=$PROJECT"-stars"
      echo "#define STARS" >> simopt.h
     ;;
   g )
      NAME=$NAME"gaussian"
      PROJECT=$PROJECT"-gaussian"
      echo "#define GAUSSIAN" >> simopt.h
     ;;
   * ) echo "ERROR: unsupported algorithm"; exit ;;
  esac

  COULOMB=0 # to ask correctly for Widom
  ;;
p )
  echo "3D periodic version"
  echo "^^^^^^^^^^^^^^^^^^^"
  echo
  . "$COOKDIR/coulomb.sh"
  ;;
s )
  echo "2D x,y-periodic (SLIT) version"
  echo "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^"
  echo
  echo "The system is NOT periodic in z."
  echo "Combination with Ewald makes it z-periodic provided that the pair interactions"
  echo "with z-neigboring periodic boxes are negligible (cutoffs < z-gap size)."
  echo

  echo "#define SLIT" >> simopt.h
  . "$COOKDIR/coulomb.sh"
  ;;
esac

echo
echo "===== Select POLARIZABILITY support:"
echo "      (polarizability requires a special version of blend and force field)"
echo "> n = No polarizability (DEFAULT)"
echo "  0 = Drude oscillators aka Charges On Springs (COS), incl. Gaussian charges,"
echo "      also scalar linear polarizability implemented via (large) point COS."
echo "------More features can be obtained by summing the following flags:"
echo "  1 = Repulsive antipolarizability (shell core), all terms included"
echo "  2 = Saturated polarizability"
echo "  4 = Use this if there are also uncharged atoms (faster code)"
echo "  8 = Axial polarizability"
# 16 = old development version: see old+misc/constrdaa.c and old+misc/gear2polg.c
echo " 32 = Fluctuating charge model (under development)"
echo " 64 = incl. intramolecular induced dipole--induced dipole interactions (ADIM)"
read
POLAR=$REPLY
if [ "$POLAR" = "" ] ; then
  POLAR="n"
fi
if [ "$POLAR" != "n" ] ; then
  echo "#define POLAR $POLAR" >> simopt.h
  NAME=$NAME"p"$POLAR
  PROJECT=$PROJECT"-polar"$POLAR
fi

echo
echo "===== Select DETAILS of boundary conditions:"
echo "> p = Plain version (DEFAULT)"
echo "  s = Slab: z-forces (incl. walls), dens.profiles, surf.tension, NPzzE melting"
echo "      also more geometries (droplet, trickle, cavities)"
echo "  m = Slab: as above and NPzzE melting"
echo "  c = Cleaving (with most Slab features)"
echo "      (see cook/generic/simopt.h for more SLAB options)"
# echo "  g = Gold (conducting) wall, charge inversion method"
# echo "      Incompatible with Ewald+linked-cell list; NOT TESTED RECENTLY"
# echo "NOTES: Versions w and g imply z-density profile measurements"
# echo "       WIDOM (if selected) is measured in z-slabs, too"
read
SV=$REPLY
if [ "$SV" = "" ] ; then
  SV="p"
fi
case "$SV" in
s )
  PROJECT=$PROJECT"-slab"
  NAME=$NAME"s"
  echo "#define SLAB 0" >> simopt.h
  ;;
m )
  PROJECT=$PROJECT"-melt"
  NAME=$NAME"s"
  echo "#define SLAB 1" >> simopt.h
  ;;
c )
  PROJECT=$PROJECT"-cleave"
  NAME=$NAME"c"
  echo "#define SLAB 2" >> simopt.h
  ;;
# w )
#   PROJECT=$PROJECT"-wall"
#   NAME=$NAME"w"
#   echo "#define WALL" >> simopt.h
#   echo "#define SLAB 0" >> simopt.h
#   ;;
# g )
#   NAME=$NAME"g"
#   PROJECT=$PROJECT"-gold"
#   echo "#define GOLD" >> simopt.h
#   echo "#define SLAB 0" >> simopt.h
#   ;;
* )
esac

echo
echo "===== The NON-BONDED forces are:"
echo "      (subdirectory of sim/ with sitesite.c and sitesite.h is specified)"
echo "      (non-lj requires special version of blend and force field)"
echo "> lj     = Lennard-Jones (DEFAULT), various combining rules available"
echo "  busing = exp-6 in the Busing/Tosi-Fumi form, C/r^6 + more"
echo "  buck   = exp-6 in A*exp(-B*r) form, C/r^6 + more (DOUBLECHECK CHANGES)"
echo "           (busing and buck are compatible with shell-core, POLAR=1)"
echo "  exp6   = exp-6 with Lennard-Jones-like combining rules, POLAR=0 supported"
echo "  metal  = tight-binding potential"
echo "  wcalj  = Lennard-Jones cut at the minimum and shifted"
echo "           NOTE: LJcutoff=-1 and corr=0 or 4 is required"
read
NONBOND=$REPLY
if [ "$NONBOND" = "" ] ; then
  NONBOND=lj
elif [ "$NONBOND" = "buck" -o "$NONBOND" = "busing" ] ; then
  echo
  echo "===== add term D/r^POW (applies to exp-6-like potentials):"
  echo "> 0 no D/r^POW added (DEFAULT)"
  echo "  8 add term -D/r^8 (Tosi-Fumi)"
  echo " 12 add term D/r^12 (repulsive core)"
  read
  POW=$REPLY
  [ "$POW" == "" ] && POW=0
  echo "#define POW $POW" >> simopt.h
fi
echo "#define NONBOND \"$NONBOND\"" >> simopt.h

if [ $BC = "f" ]; then
  MODULE=$MODULE" water"
else
  echo
  echo "===== Are all your molecules smaller than half the box size (y/N)?"
  echo "      (It is safer to say No if you are not sure)"
  read
  if [ "$REPLY" != "y" ]; then
    echo "#define WORM" >> simopt.h
  fi

  echo
  echo "===== Select algorithm for PAIR FORCES calculations (affects speed):"
  echo "      Except version \"c\", the cutoff must not exceed minimum half box."
  echo "      The cutoff refers to the elst cutoff; always cutoff>=LJcutoff."
  echo "  d = Direct evaluation of pair forces      (accepts PARALLEL=2)"
  echo "      The cutoff equals or is slightly less than half box"
  echo "  c = Direct evaluation of pair forces      (accepts PARALLEL=2)"
  echo "      The cutoff is significantly smaller than half box"
  echo "  ps = Periodic Sum                         (accepts PARALLEL=2)"
  echo "      Direct evaluation - the cutoff is longer than half box (slow!)"
  echo "> l = Linked-cell list method (DEFAULT)     (accepts PARALLEL=1 and 2)"
  echo "      Efficient for large systems and/or the cutoff sphere significantly"
  echo "      smaller than the box.  WARNING: not compatible with all other options"
  echo "  l0 = As above, with local copies of forces within cells"
  echo "  l1 = As above, with PARALLEL=1: 1st site in pair indirected, 2nd copied"
  echo "      (Which linked-cell version is more efficient depends mainly on the"
  echo "      processor and cache architecture, data size, and parallelism)"
  read
  ALG=$REPLY
  if [ "$ALG" = "" ] ; then
    ALG="l"
  fi
  case "$ALG" in
  c )
    NAME=$NAME"c"
    MODULE=$MODULE" water"
    echo "#define CUT" >> simopt.h
    ;;
  ps )
    NAME=$NAME"ps"
    MODULE=$MODULE" water"
    echo "#define PERSUM" >> simopt.h
    ;;
  l )
    MODULE=$MODULE" linklist"
    NAME=$NAME"lc"
    PROJECT=$PROJECT"-lc"
    echo "#define LINKCELL 3" >> simopt.h
    LC=3
    ;;
  l0 )
    MODULE=$MODULE" linklist"
    NAME=$NAME"lc"
    PROJECT=$PROJECT"-lc0"
    echo "#define LINKCELL 0" >> simopt.h
    LC=0
    ;;
  l1 )
    MODULE=$MODULE" linklist"
    NAME=$NAME"lc"
    PROJECT=$PROJECT"-lc1"
    echo "#define LINKCELL 1" >> simopt.h
    LC=1
    ;;
  * )
    MODULE=$MODULE" water"
    ;;
  esac

  echo
  echo "===== Select PARALLELIZATION (using linux threads):"
  echo "      (Not all combinations with the pair-forces method are valid)"
  echo "> 0 = Serial version (DEFAULT)"
  echo "  1 = Linked-cell list method parallelized in x-slabs, Ewald parallelized"
  echo " -1 = As above without barriers (deprecated: threads started/stopped more often)"
  echo "  2 = Pair forces running in parallel with Ewald (2 threads)"
  echo "      (load-balancing by hand, good for medium-size systems)"
  read
  PAR=$REPLY
  if [ "$PAR" == "0" ] ; then
# PARALLEL=0 is the same as undef, so if 0 given, no P suffix...
     PAR=""
  fi
  if [ "$PAR" != "" ] ; then
    PROJECT=$PROJECT"-par"$PAR
    NAME=$NAME"P"$PAR
    if [ "$PAR" == "-1" ] ; then
      echo "#define PARALLEL 1" >> simopt.h
    elif [ "$PAR" == "1" ] ; then
      echo "#define PARALLEL 1" >> simopt.h
      echo "#define BARRIER" >> simopt.h
    else
      echo "#define PARALLEL $PAR" >> simopt.h
    fi
  fi

  echo
  echo "===== Request PRESSURE TENSOR calculations; sum the following flags:"
  echo "  0 = no pressure tensor calculated (1-2% faster than option 3)"
  echo "  1 = calculate the virial (static) part incl. constraint forces"
  echo "  2 = calculate the kinetic part"
  echo "> 3 = 1+2 (DEFAULT)"
  echo "  4 = calculate also off-diagonal terms Pxy,Pyz,Pxz"
  echo "      (if not set, only diagonal terms Pxx,Pyy,Pzz are calculated)"
  echo "  8 = calculate also center-of-mass based (molecular) tensor"
  echo "      - applies for the kinetic part only (use with 2)!"
  echo " 16 = SPECIAL: kurtosis and 6-th normalized central moment of v_CM"
  read
  PRESSURETENSOR=$REPLY
  if [ "$PRESSURETENSOR" != "" ] ; then
    echo "#define PRESSURETENSOR $PRESSURETENSOR" >> simopt.h
    PROJECT=$PROJECT"-PT"$PRESSURETENSOR
  fi
fi

if [ "$NONBOND" = "metal" ] ; then
  echo "#define METAL" >> simopt.h
  MODULE=${MODULE/water/}
  PROJECT=$PROJECT"-METAL"
fi

echo
echo "===== Verlet algorithm (velocity, kinetic energy, pressure tensor):"
echo "  0 = v(t) = V1, Ekin = m/2 v(t)^2 (\"shifted\", in general not recommended)"
echo "  1 = v(t) = (V1+V2)/2, Ekin = m/2 v(t)^2 (Velocity Verlet)"
echo "  2 = Ekin = m/2 * V1.V2 (Ekin+Epot constant for harmonic oscillator)"
echo "> 3 = Ekin = m/4 * (V1^2+V2^2) (leap-frog, recommended default)"
echo "  5 = Ekin = m/24*(5V1^2 + 5V2^2 + 2V1.V2) (best equipartition harm.diat.)"
echo "  6 = Ekin = m/24*(7V1^2 + 7V2^2 - 2 V1/V2) (on average = Beeman)"
echo "where V1 = [r(t)-r(t-h)]/h, V2 = [r(t+h)-r(t)]/h"
echo "The formulas for the kinetic pressure tensor contain x instead of ."
read
if [ "$REPLY" = "" ]; then
  echo "#define VERLET 3" >> simopt.h
else
  echo "#define VERLET $REPLY" >> simopt.h
fi

echo "#define SHAKE 1" >> simopt.h
echo "SHAKE=1 version configured, see cook/generic/simopt.h for more options."

echo
echo "===== Select optional extensions:"
echo "      <Enter> for none, string of keys for several options"
echo "  a = anchoring a molecule and mean force calculations (ANCHOR)"
[ $BC != "f" ] && echo "  b = Bjerrum defects calculations (BJERRUM)"
echo "  c = cluster analysis (CLUSTER)"
echo "  d = dihedral (proper torsion) angle distribution (DIHEDRAL)"
echo "  f = Wang-Landau (aka metadynamics aka conformational flooding)"
[ $BC != "f" ] && echo "  e = Electronic Continuum Correction (ECC)"
echo "  h = harmonic expansion coefficients and full correlation function"
echo "      (HARMONICS) .."
[ "$LC" == "" ] && echo "  l = log separately terms as LJ, intramolecular energy (MONITOR)"
echo "  m = Morse bond potential (expanded to dist^4) support"
echo "  r = radius of gyration, end-to-end distance, and direct calculation"
echo "      of the velocity-velocity time autocorrelation functions (RGYR)"
[ $BC != "f" ] && echo "  s = periodic z-force for NEMD viscosity measurements (SHEAR)"
echo "  t = TCF of SPC water (special project)"                  
[ $COULOMB -ge "0" ] && echo "  w = Widom particle insertion method (WIDOM)"
echo "  x = geometric cross-section (XSECTION)"

read

if [ -n "$REPLY" ] ; then
  PROJECT="$PROJECT+$REPLY"
  NAME="$NAME+$REPLY"
fi

for (( i=0; i<${#REPLY}; i++ )); do
  
  case "${REPLY:$i:1}" in
    a ) echo "#define ANCHOR" >> simopt.h ;;
    b ) if [ $BC != "f" ] ; then
          echo "#define BJERRUM" >> simopt.h
          MODULE=$MODULE" minimize"
        else
          echo "BJERRUM ignored because FREE"
        fi
        ;;
    c ) echo "#define CLUSTERS" >> simopt.h ;;
    d ) echo
        echo "===== select support for DIHEDRAL (proper torsion) angle distribution measurements:"
        echo "> n = None (DEFAULT)"
        echo "  2 = Only the central bond distinguishes dihedrals"
        echo "  4 = All 4 atoms distinguish dihedrals"
        echo " -1 = All dihedrals in all molecules different:"
        echo "      Support for selection (SIMNAME.ddh) enabled"
        echo "  1 = One summary histogram for all dihedrals"
        read
        DIH=$REPLY
        if [ "$DIH" = "" ] ; then
          DIH="n"
        fi
        if [ $DIH != "n" ]; then
          echo "#define DIHHIST" $DIH >> simopt.h
        fi
        ;;
    e ) echo "#define ECC" >> simopt.h ;;
    f ) echo "#define WL" >> simopt.h ;;
    h ) echo
        echo "===== Select support for the HARMONIC EXPANSION coefficients"
        echo "      and the full correlation function:"
        echo "> n = None (DEFAULT)"
        echo "  1 = (methanol?)"
        echo "  2 = water"
        echo "  3 = linear molecule"
        read
        HARM=$REPLY
        if [ "$HARM" = "" ] ; then
          HARM="n"
        fi
        if [ $HARM != "n" ]; then
          echo "#define HARMONICS" $HARM >> simopt.h
        fi
        ;;
    l ) echo "#define LOG" >> simopt.h ;;
    m ) echo "#define MORSE" >> simopt.h ;;
    r ) echo "#define RGYR" >> simopt.h ;;
    s ) [ $BC != "f" ] && echo "#define SHEAR" >> simopt.h || echo "SHEAR ignored because FREE" ;;
    t ) echo "#define SPCTCF" >> simopt.h ;; 
    w ) [ $COULOMB -ge "0" ] && echo "#define WIDOM" >> simopt.h || echo "WIDOM ignored because COULOMB=$COULOMB" ;;
    x ) echo "#define XSECTION" >> simopt.h ;;
  esac
done

echo "#define PROJECT" \"$PROJECT\" >> simopt.h

if [ -e $METAMAKE ] ; then
  mv $METAMAKE $METAMAKE~
  echo "NOTE: existing (previous) $METAMAKE was renamed to $METAMAKE~"
  echo
fi

echo "!! makemake $CC" > $METAMAKE
echo "!! This file was generated by configure.sh" >> $METAMAKE
echo "!! See cook/generic/$METAMAKE for more details" >> $METAMAKE
echo >> $METAMAKE
echo "O := -DSCR" >> $METAMAKE

if [ -e ../home.mmk ]; then
  echo "!include \"../home.mmk\"" >> $METAMAKE
elif [ -e ../../home.mmk ]; then
  echo "!include \"../../home.mmk\"" >> $METAMAKE
elif [ -e ../../../home.mmk ]; then
  echo "!include \"../../../home.mmk\"" >> $METAMAKE
elif [ -e ../../../../home.mmk ]; then
  echo "!include \"../../../../home.mmk\"" >> $METAMAKE
else
  echo "ERROR: home.mmk not found!"
  echo "(none of ../../home.mmk ../../../home.mmk ../../../../home.mmk exists)"
  echo "The directory in which cook is configured must be the 1st, 2nd"
  echo "or 3rd-level subdirectry of MACSIMUS/cook/."
  exit
fi

echo "!! project directories, relative to !home (see \"home.mmk\")" >> $METAMAKE
echo "!dir = gen sim sim/$NONBOND cook cook/$SUBDIR" >> $METAMAKE

if [ "$EXENAME" != "" ] ; then
  NAME="$EXENAME"
else
  echo "The executable name was determined from the configuration"
fi

echo "!! project dependencies, will be expanded to makefile by \"makemake\"" >> $METAMAKE
echo -n "!"$NAME ": " >> $METAMAKE
echo "ground statics varfile bitfile cputime conjgrad gjlineq asksig \\" >> $METAMAKE
echo "  simcg simgear interpot setss intrapot sitesite norm xforces linregr \\" >> $METAMAKE
echo "  constrd rhs simglob pakcp simils siminit main simpot simmeas simmeasx \\" >> $METAMAKE
echo "  forces simdef \\" >> $METAMAKE
echo " "$MODULE >> $METAMAKE
echo >> $METAMAKE

echo "The executable name is \"$PWD/$NAME\""
echo "You may change it in file \"$PWD/$METAMAKE\""
echo
echo "<Enter> to convert $METAMAKE to the makefile (linux version), ^C to cancel"

read

if [ -e ../../bin/makemake ]; then
  MM="../../bin/makemake"
else
  MM=makemake
fi

if [ "$PAR" != "" ] ; then
  MMOPT=pthread
fi

echo $MM $CC $MMOPT
if ! $MM $CC $MMOPT ; then
  echo $MM $CC $MMOPT failed
  exit 1
fi

echo
echo -n "SITE-SITE: "
fgrep "D3 =" makefile

echo "<Enter> to erase all previous *.o modules and compile, ^C to cancel"
read
rm *.o

echo "========================================================================"
make "$NAME"
echo "recompile by:"
echo "  cd $SUBDIR; make $NAME"
