#!/usr/bin/env bash
# (note: Classic Bourne shell (#!/bin/sh) also worked at some point.)

# Author: Andrew Jewett (jewett.aij at g mail)
# License: MIT License  (See LICENSE.md)
# Copyright (c) 2013

G_PROGRAM_NAME="moltemplate.sh"
G_VERSION="2.22.4"
G_DATE="2025-8-21"

echo "${G_PROGRAM_NAME} v${G_VERSION} ${G_DATE}" >&2
echo "" >&2

# Check for python:

if which python3 > /dev/null; then
    PYTHON_COMMAND='python3'
elif which python > /dev/null; then
    PYTHON_COMMAND='python'
    if [ `python --version | awk '{print substr($2,0,2)}'` -lt 3 ]; then
        echo "############################################################" >&2
        echo "############################################################" >&2
        echo "############################################################" >&2
        echo "  WARNING: Support for python2.7 has been depreciated." >&2
        echo "           Please install python version 3.4 or later." >&2
        echo "           Use older versions of python at your own risk." >&2
        echo "############################################################" >&2
        echo "############################################################" >&2
        echo "############################################################" >&2
    fi
elif which python2 > /dev/null; then
    PYTHON_COMMAND='python2'
    echo "############################################################" >&2
    echo "############################################################" >&2
    echo "############################################################" >&2
    echo "  WARNING: Support for python2.7 has been depreciated." >&2
    echo "           Please install python version 3.4 or later." >&2
    echo "           Use older versions of python at your own risk." >&2
    echo "############################################################" >&2
    echo "############################################################" >&2
    echo "############################################################" >&2
else
    echo "Error:  $G_PROGRAM_NAME requires python, python3, or python2." >&2
    exit 1
fi


# First, determine the directory in which the python scripts are located.
# (such as ttree.py).  It could either be the directory where the script
# file is located, OR it could be the parent of this directory.
PY_SCR_DIR=`dirname "$0"`
if [ ! -s "${PY_SCR_DIR}/ttree.py" ]; then
    PY_SCR_DIR="$PY_SCR_DIR/.."
fi


MSG_BAD_INSTALL=$(cat <<EOF

INSTALLATION ERROR:
Follow the instructions in the "Installation" chapter of the moltemplate manual.
(Note: You may need to log out and log in again before the changes take effect.)

EOF
)


# Maintain a list of papers that we would like users to cite.
# (This list currently has only one entry.)
export MOLTEMPLATE_CITE_LIST="Jewett et al. J.Mol.Biol. (2021) (https://doi.org/10.1016/j.jmb.2021.166841)"



ERR_BAD_INSTALL()
{
  echo "$MSG_BAD_INSTALL" >&2
  exit 1
}


ERR_INTERNAL()
{
  echo "    !!!!!!   Possible internal error  !!!!!!" >&2
  echo "This could be a bug in moltemplate." >&2
  echo "Please report this error." >&2
  echo "(And please include the last few lines of moltemplate output preceeding this.)" >&2
  echo "  Thank you." >&2
  exit 100
}


MOLTEMPLATE_FILES_NEEDED=$(cat <<EOF
ttree.py
lttree.py
lttree_check.py
lttree_postprocess.py
nbody_by_type.py
nbody_fix_ttree_assignments.py
nbody_reorder_atoms.py
pdbsort.py
postprocess_input_script.py
postprocess_coeffs.py
remove_duplicate_atoms.py
remove_duplicates_nbody.py
renumber_DATA_first_column.py
ttree_render.py
dump2data.py
raw2data.py
EOF
)

# Save the default Internal Field Separator and define "newline" as an alternative IFS.
OIFS=$IFS
#CR=$'\n'
CR="
"

IFS=$CR
for f in $MOLTEMPLATE_FILES_NEEDED; do
    if [ ! -s "${PY_SCR_DIR}/$f" ]; then
        echo "Error: Missing file \"${PY_SCR_DIR}/$f\"" >&2
        ERR_BAD_INSTALL
    fi
done
IFS=$OIFS

# command that invokes lttree.py
LTTREE_COMMAND="$PYTHON_COMMAND \"${PY_SCR_DIR}/lttree.py\""

# command that invokes lttree_check.py
LTTREE_CHECK_COMMAND="$PYTHON_COMMAND \"${PY_SCR_DIR}/lttree_check.py\""

# command that invokes lttree_postprocess.py
LTTREE_POSTPROCESS_COMMAND="$PYTHON_COMMAND \"${PY_SCR_DIR}/lttree_postprocess.py\""


# -----------------------------------------------------------
# If everything worked, then running ttree usually
# generates the following files:
#
# Users of lttree typically generate the following files:
# The variable below refer to file names generated by
# write() and write_once() commands in a lttree-file.
# (I keep changing my mind what I want these names to be.)
data_prefix="Data "
data_prefix_no_space="Data"
data_header="Data Header"
data_atoms="Data Atoms"
data_masses="Data Masses"
data_velocities="Data Velocities"
data_bonds="Data Bonds"
data_bond_list="Data Bond List"
data_bonds_atomid_atomid="Data Bonds AtomId AtomId"
data_angles="Data Angles"
data_dihedrals="Data Dihedrals"
data_impropers="Data Impropers"
data_bond_coeffs="Data Bond Coeffs"
data_angle_coeffs="Data Angle Coeffs"
data_dihedral_coeffs="Data Dihedral Coeffs"
data_improper_coeffs="Data Improper Coeffs"
data_pair_coeffs="Data Pair Coeffs"
data_pairij_coeffs="Data PairIJ Coeffs"
data_cmap="Data CMAP"

# interactions-by-type (not id. This is not part of the LAMMPS standard.)
data_charge_by_bond="Data Charge By Bond"
data_bonds_by_type="Data Bonds By Type"
data_angles_by_type="Data Angles By Type"
data_dihedrals_by_type="Data Dihedrals By Type"
data_impropers_by_type="Data Impropers By Type"

# class2 data sections
data_bondbond_coeffs="Data BondBond Coeffs"
data_bondangle_coeffs="Data BondAngle Coeffs"
data_middlebondtorsion_coeffs="Data MiddleBondTorsion Coeffs"
data_endbondtorsion_coeffs="Data EndBondTorsion Coeffs"
data_angletorsion_coeffs="Data AngleTorsion Coeffs"
data_angleangletorsion_coeffs="Data AngleAngleTorsion Coeffs"
data_bondbond13_coeffs="Data BondBond13 Coeffs"
data_angleangle_coeffs="Data AngleAngle Coeffs"

# sections for non-point-like particles:
data_ellipsoids="Data Ellipsoids"
data_lines="Data Lines"
data_triangles="Data Triangles"

# periodic boundary conditions
data_boundary="Data Boundary"
# (for backward compatibility), an older version of this file was named:
data_pbc="Data PBC"

# ---------------------------------------------------------------
# Note: The files above are fragments of a LAMMPS data file.
# In addition, moltemplate will probably also generate the following files:
# (These files represent different sections of the LAMMPS input script.)
# ---------------------------------------------------------------

in_prefix="In "
in_prefix_no_space="In"
in_init="In Init"
in_settings="In Settings"
in_coords="In Coords"
in_charges="In Charges"

#     If present, the various "In " files contain commands which should be
#     included by the user in their LAMMPS input script. (This task is left
#     to the user.) However, the "Data " files are processed and pasted together
#     automatically in order to build a LAMMPS data file.
# ---------------------------------------------------------------

tmp_atom_coords="tmp_atom_coords.dat"  #<-temporary file for storing coordinates
tmp_dump="tmp_dump.dat"
tmp_ellips_quat="tmp_ellips_quat.dat"


MOLTEMPLATE_TEMP_FILES=$(cat <<EOF
*.template
ttree_assignments.txt
$tmp_atom_coords
$tmp_dump
$tmp_ellips_quat
$data_masses
$data_pair_coeffs
$data_pairij_coeffs
$data_bond_coeffs
$data_angle_coeffs
$data_dihedral_coeffs
$data_improper_coeffs
$data_atoms
$data_velocities
$data_bonds
$data_bond_list
$data_angles
$data_dihedrals
$data_impropers
$data_bondbond_coeffs
$data_bondangle_coeffs
$data_middlebondtorsion_coeffs
$data_endbondtorsion_coeffs
$data_angletorsion_coeffs
$data_angleangletorsion_coeffs
$data_bondbond13_coeffs
$data_angleangle_coeffs
$data_ellipsoids
$data_lines
$data_triangles
$data_boundary
$data_header
$data_bonds_by_type*
${data_angles_by_type}*
${data_dihedrals_by_type}*
${data_impropers_by_type}*
$data_charge_by_bond
$in_init
$in_settings
EOF
)


IFS=$CR
for f in $MOLTEMPLATE_TEMP_FILES; do
    #echo "removing [$f]"
    rm -f "$f"
done
IFS=$OIFS

rm -rf output_ttree




SYNTAX_MSG=$(cat <<EOF
Syntax example:
Usage:

moltemplate.sh [-atomstyle style] \
               [-pdb/-xyz/-dump/-raw coord_file] \
               [-a assignments.txt] file.lt

Optional arguments:

-atomstyle style  By default, moltemplate.sh assumes you are using the "full"
                atom style in LAMMPS.  You can change the atom style to "dipole"
                using -atomstyle dipole.  If you are using a hybrid style,
                you must enclose the list of styles in quotes.  For example:
                -atomstyle "hybrid full dipole"
                For custom atom styles, you can also specify the
                list of column names manually (enclosed in quotes):
                -atomstyle "molid x y z atomid atomtype mux muy muz"

-xyz xyz_file   An optional xyz_file argument can be supplied as an argument
                following "-xyz".
                This file should contain the atomic coordinates in xyz format.
                (The atoms must appear in the same order in the data file.)

-pdb pdb_file   An optional pdb_file argument can be supplied as an argument
                following "-pdb".

                This should be a PDB file (with ATOM or HETATM records) with
                the coordinates you wish to appear in the LAMMPS data file.
                (The atoms must appear in the same order in the data file.)

                If the PDB file contains periodic boundary box information
                (ie., a "CRYST1" record), this information is also copied
                to the LAMMPS data file.
                (Support for triclinic cells is experimental as of 2012-2-13.
                 Other molecular structure formats may be supported later.)

-dump dump_file  An optional dump_file argument can be supplied as an argument
                following "-dump".

                This option should be used to overwrite the coordinate section
                in the LAMMPS data file. If present, it will also write the
                Ellipsoids and Velocities sections.

                For data style ellipsoid, the dump file MUST include the
                quaterions printed using one of the following computes:
                compute q      all property/atom quatw quati quatj quatk
		compute orient all property/atom quati quatj quatk quatw
		The parser also accepts DUMP headers with "edited" names
		such as qi, qj, qk, qw and quati, quatj, quatk, quatw.

                The box boundaries in the DUMP are used in the LAMMPS data file.

-raw raw_file   An optional file containing RAW atomic coordinates. RAW files are
                simple 3-column ASCII files contain X Y Z coordinates for every
		atom, separated by spaces.
		

-a "@atom:x 1"
-a assignments.txt
                The user can customize the numbers assigned to atom, bond,
                angle, dihedral, and improper types or id numbers by using
                   -a "VARIABLE_NAME VALUE"
                for each variable you want to modify.  If there are many
                variables you want to modify, you can save them in a file
                (one variable per line).  For an example of the file format
                run moltemplat.sh once and search for a file named
                "ttree_assignments.txt".  (This file is often located in
                the "output_ttree/" directory.) Once assigned, the remaining
                variables in the same category will be automatically assigned
                to values which do not overlap with your chosen values.
-b assignments.txt
                "-b" is similar to "-a". However, in this case, no attempt
                is made to assign exclusive (unique) values to each variable.
-nocheck
                Normally moltemplate.sh checks for common errors and typos and
                halts if it thinks it has found one.  This forces the variables
                and categories as well as write(file) and write_once(file)
                commands to obey standard naming conventions.  The "-nocheck"
                argument bypasses these checks and eliminates these restrictions.

-checkff        This cause moltemplate.sh to check to make sure that there
                are valid angle and dihedral interactions defined for every
                3 or 4 consecutively bonded atoms in the system
                (defined in "Angles/Dihedrals By Type").

-overlay-bonds      Normally, moltemplate.sh checks to see if multiple angle
-overlay-angles     interactions are defined for the same triplet of atoms.
-overlay-dihedrals  If so, it deletes the redundant ones (keeping the last one).
-overlay-impropers  (It does the same thing for bonds, dihedrals, and impropers.)
-overlay-all        Use these options to prevent that behavoir, or to prevent
                    reordering of atom indexes, e.g. by allowing to define a bond
                    between atom 3 and 2 instead of 2 and 3.

-report-duplicates filter_id filter_type
              Warn when multiple bonds, angles, dihedrals, or impropers were found
              between the same atoms (regardless of whether ``-overlay-'' was used).
              The filter_id and filter_type arguments will ignore interactions whose
              ids or types lack those strings. (Normally set them to "".  See manual.)

-angle-symmetry file.py     Normally moltemplate.sh reorders the atoms in each
-dihedral-symmetry file.py  angle, dihedral, improper, and bond interaction.
-improper-symmetry file.py  TO TURN OFF ATOM REORDERING, SET file.py to "NONE"
-bond-symmetry file.py      You can override the rules that moltemplate.sh uses
                            by supplying a file (file.py) with symmetry rules.
                            See nbody_Dihedrals.py, nbody_Impropers.py (in the
                            moltemplate directory) to learn the file format.

-molc           Additional post-processing for the file "In Settings". This
                options implicitly sets "-overlay-bonds".

-lammps-script basename 
                By default, moltemplate.sh writes a LAMMPS execution script with
		the same base name of the input LT file. If the LT file has the
		special name "system.lt", then the LAMMPS execution script is
		named run.in.EXAMPLE. This behaviour can be overwritten by
		specifying a custom basename for the execution script. Note that
		the INIT, DATA, SETTINGS, and RUN files will still have the same
		basename of the parent LT file.

-run-example    Append an example run section at the end of the example input
                script. The commands are commented out and need manual editing
		to get a meaningful simulation.

-forbid-wildcards  Forbid the use of "*" and "?" characters in "pair coeff",
                "bond coeff", "angle coeff", "dihedral coeff", and
		"improper coeff" commands.

EOF
)

# --- Periodic boundary box information (default) ---

# We will determine these numbers later.
TRICLINIC=""
BOXSIZE_MINX=0.0
BOXSIZE_MINY=0.0
BOXSIZE_MINZ=0.0
BOXSIZE_MAXX=""
BOXSIZE_MAXY=""
BOXSIZE_MAXZ=""
BOXSIZE_XY=0.0
BOXSIZE_XZ=0.0
BOXSIZE_YZ=0.0



if [ "$1" = "--help" ]; then
    echo "$SYNTAX_MSG" >&2
    exit 0
fi



# --- Did the user specify a file containing atomic coordinates?

rm -f "$tmp_atom_coords" "$tmp_dump" "$tmp_ellips_quat"

# Optional files containing atom coordinates:
PDB_FILE=""
XYZ_FILE=""
DUMP_FILE=""
RAW_FILE=""
LT_FILE=""
OUT_FILE_BASE="system"
OUT_FILE_EXAMPLE_SCRIPT=""
# REMOVE_DUPLICATE variables:
# ...If true (default), then any duplicate entries in the lists of bonds
# bonds, angles, dihedrals, or impropers in the LAMMPS DATA file
# are removed, giving priority to the most recent entry in the list.
# (This might not be necessary, but I want to be careful.)
REMOVE_DUPLICATE_BONDS="true"
REMOVE_DUPLICATE_ANGLES="true"
REMOVE_DUPLICATE_DIHEDRALS="true"
REMOVE_DUPLICATE_IMPROPERS="true"
REPORT_DUPLICATES=""
RD_ID_FILTER=""
RD_TYPE_FILTER=""
SETTINGS_MOLC=""
CHECKFF=""
RUN_VMD_AT_END=""
APPEND_EXAMPLE_SCRIPT=""


ARGC=0
for A in "$@"; do
    A_FIRSTCHAR="$(echo $A| cut -c 1)"
    # (Note to self: this next line only works in bash, not classic sh.)
    if [ "$A_FIRSTCHAR" = "\$" ]; then
        A="\\$A" # put an extra slash in front to prevent expansion later
    fi
    ARGC=$((ARGC+1))
    eval ARGV${ARGC}=\"$A\"
done

TTREE_ARGS=""
ATOM_STYLE="full"
ATOM_STYLE_ARG="-atomstyle \"$ATOM_STYLE\""

i=0
while [ "$i" -lt "$ARGC" ]; do
    i=$((i+1))
    eval A=\${ARGV${i}}

    if [ "$A" = "-nocheck" ]; then
        # Disable syntax checking by undefining LTTREE_CHECK_COMMAND
        unset LTTREE_CHECK_COMMAND
        unset LTTREE_POSTPROCESS_COMMAND
    elif [ "$A" = "-allow-wildcards" ]; then
        if [ -z "$LTTREE_CHECK_ARGS" ]; then
            LTTREE_CHECK_ARGS="\"$A\""
        else
            LTTREE_CHECK_ARGS="${LTTREE_CHECK_ARGS} \"$A\""
        fi
    elif [ "$A" = "-forbid-wildcards" ]; then
        if [ -z "$LTTREE_CHECK_ARGS" ]; then
            LTTREE_CHECK_ARGS="\"$A\""
        else
            LTTREE_CHECK_ARGS="${LTTREE_CHECK_ARGS} \"$A\""
        fi
    elif [ "$A" = "-checkff" ]; then
        CHECKFF="$A"
    elif [ "$A" = "-overlay-bonds" ]; then
        # In that case, do not remove duplicate bond interactions
        unset REMOVE_DUPLICATE_BONDS
    elif [ "$A" = "-overlay-angles" ]; then
        # In that case, do not remove duplicate angle interactions
        unset REMOVE_DUPLICATE_ANGLES
    elif [ "$A" = "-overlay-dihedrals" ]; then
        # In that case, do not remove duplicate dihedral interactions
        unset REMOVE_DUPLICATE_DIHEDRALS
    elif [ "$A" = "-overlay-impropers" ]; then
        # In that case, do not remove duplicate improper interactions
        unset REMOVE_DUPLICATE_IMPROPERS
    elif [ "$A" = "-overlay-all" ]; then
        # In that case, do not remove duplicate bonded interactions, nor reorder the bonded terms.
        unset REMOVE_DUPLICATE_BONDS
        unset REMOVE_DUPLICATE_ANGLES
        unset REMOVE_DUPLICATE_DIHEDRALS
        unset REMOVE_DUPLICATE_IMPROPERS
    elif [ "$A" = "-report-duplicates" ]; then
        REPORT_DUPLICATES="true"
        if [ "$((i+1))" -eq "$ARGC" ]; then
            echo "ERROR: Expected 2 string arguments following the -report-duplicates argument" >&2
            exit 7
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        RD_ID_FILTER="$A"
        i=$((i+1))
        eval A=\${ARGV${i}}
        RD_TYPE_FILTER="$A"
    elif [ "$A" = "-vmd" ]; then
        RUN_VMD_AT_END="true"
    elif [ "$A" = "-molc" ]; then
        # Set the -overlay-bonds, if not specified otherwise.
        unset REMOVE_DUPLICATE_BONDS
        SETTINGS_MOLC="true"
    elif [ "$A" = "-run-example" ]; then
        APPEND_EXAMPLE_SCRIPT="true"
    elif [ "$A" = "-lammps-script" ]; then
        if [ "$i" -eq "$ARGC" ]; then
            echo "ERROR: File name expected following the -lammps-script argument" >&2
            exit 7
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        OUT_FILE_EXAMPLE_SCRIPT="$A"
    elif [ "$A" = "-raw" ]; then
        if [ "$i" -eq "$ARGC" ]; then
            echo "$SYNTAX_MSG" >&2
            exit 7
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        RAW_FILE=$A
        if [ ! -s "$RAW_FILE" ]; then
            echo "$SYNTAX_MSG" >&2
            echo "-----------------------" >&2
            echo "" >&2
            echo "Error: Unable to open RAW-file \"$RAW_FILE\"." >&2
            echo "       (File is empty or does not exist.)" >&2
            exit 8
        fi
        #echo "  (extracting coordinates from \"$RAW_FILE\")" >&2
        awk '{if (NF==3) {print $0}}' < "$RAW_FILE" > "$tmp_atom_coords"

    elif [ "$A" = "-bond-symmetry" ]; then
        # Change the atom ordering rules in a 2-body bonded interaction:
        if [ "$i" -eq "$ARGC" ]; then
            echo "$SYNTAX_MSG" >&2
            exit 7
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        if [ "$A" = "NONE" ] || [ "$A" = "none" ] || [ "$A" = "None" ]; then
            SUBGRAPH_SCRIPT_BONDS="bonds_nosym.py"
        else
            SUBGRAPH_SCRIPT_BONDS=$A
        fi

    elif [ "$A" = "-angle-symmetry" ]; then
        # Change the atom ordering rules in a 3-body "angle" interaction:
        if [ "$i" -eq "$ARGC" ]; then
            echo "$SYNTAX_MSG" >&2
            exit 7
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        if [ "$A" = "NONE" ] || [ "$A" = "none" ] || [ "$A" = "None" ]; then
            SUBGRAPH_SCRIPT_ANGLES="angles_nosym.py"
        else
            SUBGRAPH_SCRIPT_ANGLES=$A
        fi

    elif [ "$A" = "-dihedral-symmetry" ]; then
        # Change the atom ordering rules in a 4-body "dihedral" interaction:
        if [ "$i" -eq "$ARGC" ]; then
            echo "$SYNTAX_MSG" >&2
            exit 7
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        if [ "$A" = "NONE" ] || [ "$A" = "none" ] || [ "$A" = "None" ]; then
            SUBGRAPH_SCRIPT_DIHEDRALS="dihedrals_nosym.py"
        else
            SUBGRAPH_SCRIPT_DIHEDRALS=$A
        fi

    elif [ "$A" = "-improper-symmetry" ]; then
        # Change the atom ordering rules in a 4-body "improper" interaction:
        if [ "$i" -eq "$ARGC" ]; then
            echo "$SYNTAX_MSG" >&2
            exit 7
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        if [ "$A" = "NONE" ] || [ "$A" = "none" ] || [ "$A" = "None" ]; then
            SUBGRAPH_SCRIPT_IMPROPERS="impropers_nosym.py"
        else
            SUBGRAPH_SCRIPT_IMPROPERS=$A
        fi

    elif [ "$A" = "-xyz" ]; then
        if [ "$i" -eq "$ARGC" ]; then
            echo "$SYNTAX_MSG" >&2
            exit 7
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        XYZ_FILE=$A
        if [ ! -s "$XYZ_FILE" ]; then
            echo "$SYNTAX_MSG" >&2
            echo "-----------------------" >&2
            echo "" >&2
            echo "Error: Unable to open XYZ-file \"$XYZ_FILE\"." >&2
            echo "       (File is empty or does not exist.)" >&2
            exit 8
        fi
        #echo "  (extracting coordinates from \"$XYZ_FILE\")" >&2

        # "isnum(x)" returns 1 or 0 depending upon whether or not
        # a string is numeric.
        #http://rosettacode.org/wiki/Determine_if_a_string_is_numeric#AWK

        awk '{if (NR>2) {if (NF==4) {print $2" "$3" "$4} else {if (NF==3) {print $1" "$2" "$3}}}}' < "$XYZ_FILE" > "$tmp_atom_coords"

    elif [ "$A" = "-pdb" ]; then
        if [ "$i" -eq "$ARGC" ]; then
            echo "$SYNTAX_MSG" >&2
            exit 9
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        PDB_FILE=$A
        if [ ! -s "$PDB_FILE" ]; then
            echo "$SYNTAX_MSG" >&2
            echo "-----------------------" >&2
            echo "" >&2
            echo "Error: Unable to open PDB-file \"$PDB_FILE\"." >&2
            echo "       (File is empty or does not exist.)" >&2
            exit 10
        fi
        #echo "  (extracting coordinates from \"$PDB_FILE\")" >&2
        if grep -q '^ATOM  \|^HETATM' "$PDB_FILE"; then
            # Extract the coordinates from the PDB file:

            # COMMENT:
            # I used to sort the PDB file by (ChainID,SeqNum,InsertCode)
            # and then extract the coordinates from the file.
            # This turned out to be inconvenient for users.  Instead
            # just read the coordinates in the order they appear in the file:

            awk '/^ATOM  |^HETATM/{print substr($0,31,8)" "substr($0,39,8)" "substr($0,47,8)}' \
                < "$PDB_FILE" \
                > "$tmp_atom_coords"

        else
            echo "$SYNTAX_MSG" >&2
            echo "-----------------------" >&2
            echo "" >&2
            echo "Error: File \"$PDB_FILE\" is not a valid PDB file." >&2
            exit 11
        fi
        # Now extract the periodic bounding-box informatio from the PDB file
        # The CRYST1 records are described at:
        # http://deposit.rcsb.org/adit/docs/pdb_atom_format.html
        BOXSIZE_A=-1.0
        BOXSIZE_B=-1.0
        BOXSIZE_C=-1.0
        ALPHA=" 90.00"  #Note: The space before the number in " 90.00" is intentional.
        BETA=" 90.00"   #      Later we will check to see if the system is triclinic
        GAMMA=" 90.00"  #      by comparing these strings for equality with " 90.00"
        if grep -qF "CRYST1" "$PDB_FILE"; then
            BOXSIZE_A=`awk '/CRYST1/{print substr($0,8,8)}' < "$PDB_FILE"`
            BOXSIZE_B=`awk '/CRYST1/{print substr($0,17,8)}' < "$PDB_FILE"`
            BOXSIZE_C=`awk '/CRYST1/{print substr($0,26,8)}' < "$PDB_FILE"`
            ALPHA=`awk '/CRYST1/{print substr($0,35,6)}' < "$PDB_FILE"`
            BETA=`awk '/CRYST1/{print substr($0,42,6)}' < "$PDB_FILE"`
            GAMMA=`awk '/CRYST1/{print substr($0,49,6)}' < "$PDB_FILE"`
        fi
        if [ `echo "$ALPHA!=90.0"|bc` -eq 1 ] || [ `echo "$BETA!=90.0"|bc` -eq 1 ] || [ `echo "$GAMMA!=90.0"|bc` -eq 1 ]; then
            # I transform the parameters from one format to the other by inverting
            # the transformation formula from the LAMMPS documentation 
            # https://docs.lammps.org/Howto_triclinic.html  (which matches
            # http://www.ccl.net/cca/documents/molecular-modeling/node4.html)
            TRICLINIC="true"
            PI=3.1415926535897931
            BOXSIZE_X=$BOXSIZE_A
            BOXSIZE_Y=`awk -v PI="$PI" -v BOXSIZE_B="$BOXSIZE_B" -v GAMMA="$GAMMA" 'BEGIN{print BOXSIZE_B*sin(GAMMA*PI/180.0)}'`
            BOXSIZE_XY=`awk -v PI="$PI" -v BOXSIZE_B="$BOXSIZE_B" -v GAMMA="$GAMMA" 'BEGIN{print BOXSIZE_B*cos(GAMMA*PI/180.0)}'`
            BOXSIZE_XZ=`awk -v PI="$PI" -v BOXSIZE_C="$BOXSIZE_C" -v BETA="$BETA"   'BEGIN{print BOXSIZE_C*cos(BETA*PI/180.0)}'`
            BOXSIZE_YZ=`awk -v PI="$PI" -v BOXSIZE_C="$BOXSIZE_C" -v ALPHA="$ALPHA" -v BETA="$BETA" -v GAMMA="$GAMMA" 'BEGIN{ca=cos(ALPHA*PI/180.0); cb=cos(BETA*PI/180.0); cg=cos(GAMMA*PI/180.0); sg=sin(GAMMA*PI/180.0); print BOXSIZE_C*(ca-(cg*cb))/sg}'`
            BOXSIZE_Z=`awk -v PI="$PI" -v BOXSIZE_C="$BOXSIZE_C" -v ALPHA="$ALPHA" -v BETA="$BETA" -v GAMMA="$GAMMA" 'BEGIN{ca=cos(ALPHA*PI/180.0); cb=cos(BETA*PI/180.0); cg=cos(GAMMA*PI/180.0); sg=sin(GAMMA*PI/180.0); print BOXSIZE_C*sqrt(1.0+2*ca*cb*cg-ca*ca-cb*cb-cg*cg)/sg}'`
            #BOXSIZE_Z=`awk -v BOXSIZE_C="$BOXSIZE_C" -v BOXSIZE_XZ="$BOXSIZE_XZ" -v BOXSIZE_YZ="$BOXSIZE_YZ" 'BEGIN{print sqrt((BOXSIZE_C*BOXSIZE_C)-((BOXSIZE_XZ*BOXSIZE_XZ)+(BOXSIZE_YZ*BOXSIZE_YZ)))}'`
            # Note: The two expressions for "BOXSIZE_Z" are algebraically equivalent.
        else
            BOXSIZE_X=$BOXSIZE_A
            BOXSIZE_Y=$BOXSIZE_B
            BOXSIZE_Z=$BOXSIZE_C
            BOXSIZE_XY=0.0
            BOXSIZE_XZ=0.0
            BOXSIZE_YZ=0.0
        fi
        BOXSIZE_MINX=0.0
        BOXSIZE_MINY=0.0
        BOXSIZE_MINZ=0.0
        BOXSIZE_MAXX=$BOXSIZE_X
        BOXSIZE_MAXY=$BOXSIZE_Y
        BOXSIZE_MAXZ=$BOXSIZE_Z

    # Contributing author for read DUMP: Otello M Roscioni.
    elif [ "$A" = "-dump" ]; then
        if [ "$i" -eq "$ARGC" ]; then
            echo "$SYNTAX_MSG" >&2
            exit 7
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        DUMP_FILE=$A
        if [ ! -s "$DUMP_FILE" ]; then
            echo "$SYNTAX_MSG" >&2
            echo "-----------------------" >&2
            echo "" >&2
            echo "Error: Unable to open DUMP-file \"$DUMP_FILE\"." >&2
            echo "       (File is empty or does not exist.)" >&2
            exit 8
        fi

        # Reorder the last frame of a DUMP file, as the atoms are dumped randomly.
        last=$(head -4 "$DUMP_FILE" | awk 'NR==4{print $1+9}')
        tail -$last "$DUMP_FILE" | head -9 > "$tmp_dump"
        tail -$last "$DUMP_FILE" | tail -n+10 | sort -n >> "$tmp_dump"

        # Read the box
        IFS=$CR
        box=( $(head -8 "$tmp_dump" | awk 'NR>5{for (i=1;i<=NF;i++){printf "%g\n",$i}}') )

        # Find the columns of: position, quaternion, velocity, and angular momentum.
        pos=( $(sed -n '9p;9q' "$tmp_dump" | awk '{for(i=1; i<=NF; i++){if($i~/^[xyz]/){printf "%i\n",i-2} }}') )
        # Quaternion order:
	# (I,J,K,W) from "compute orient all property/atom quati quatj quatk quatw"
	# (W,I,J,K) from "compute q      all property/atom quatw quati quatj quatk"
	quat=( $(sed -n '9p;9q' "$tmp_dump" | awk '{for(i=1; i<=NF; i++){
               if(                         $i~/^qw$/||$i~/^quatw$/||$i~/^c_q\[1\]$/||$i~/^c_orient\[4\]$/){qo[1]=i-2}
               if($i~/^qi$/||$i~/^quati$/||$i~/^qx$/||$i~/^quatx$/||$i~/^c_q\[2\]$/||$i~/^c_orient\[1\]$/){qo[2]=i-2}
               if($i~/^qj$/||$i~/^quatj$/||$i~/^qy$/||$i~/^quaty$/||$i~/^c_q\[3\]$/||$i~/^c_orient\[2\]$/){qo[3]=i-2}
               if($i~/^qk$/||$i~/^quatk$/||$i~/^qz$/||$i~/^quatz$/||$i~/^c_q\[4\]$/||$i~/^c_orient\[3\]$/){qo[4]=i-2}
               }}END{for(i in qo){printf "%i\n",qo[i]}}') )
        vel=( $(sed -n '9p;9q' "$tmp_dump" | awk '{for(i=1; i<=NF; i++){if($i~/v[xyz]/){printf "%i\n",i-2} }}') )
        angmom=( $(sed -n '9p;9q' "$tmp_dump" | awk '{for(i=1; i<=NF; i++){if($i~/angmom[xyz]/||$i~/AngularMomentum[XYZ]/){printf "%i\n",i-2} }}') )
        IFS=$OIFS

        # Orthorombic box.
        if [ ${#box[@]} == 6 ]; then
          BOXSIZE_MINX=${box[0]}
          BOXSIZE_MAXX=${box[1]}
          BOXSIZE_MINY=${box[2]}
          BOXSIZE_MAXY=${box[3]}
          BOXSIZE_MINZ=${box[4]}
          BOXSIZE_MAXZ=${box[5]}
        fi

        # Triclinic box.
        # For triclinic systems the first two columns of "ITEM: BOX BOUNDS" in 
        # the dump describe a bounding box around the system, not the system 
        # itself. See https://docs.lammps.org/dump.html and https://docs.lammps.org/Howto_triclinic.html
        if [ ${#box[@]} == 9 ]; then
          xtilt=( $(echo "0.0 ${box[2]} ${box[5]} $(bc -l<<<${box[2]}+${box[5]})"|xargs -n1|sort -n) )
          ytilt=( $(echo "0.0 ${box[8]}"|xargs -n1|sort -n) )
          BOXSIZE_MINX=$(echo ${box[0]}" "${xtilt[0]} |awk '{print $1-$2}')
          BOXSIZE_MAXX=$(echo ${box[1]}" "${xtilt[3]} |awk '{print $1-$2}')
          BOXSIZE_MINY=$(echo ${box[3]}" "${ytilt[0]} |awk '{print $1-$2}')
          BOXSIZE_MAXY=$(echo ${box[4]}" "${ytilt[1]} |awk '{print $1-$2}')
          BOXSIZE_MINZ=${box[6]}
          BOXSIZE_MAXZ=${box[7]}
          BOXSIZE_XY=${box[2]}
          BOXSIZE_XZ=${box[5]}
          BOXSIZE_YZ=${box[8]}
          TRICLINIC="true"
        fi

        # Save the coordinates.
        awk -v x=${pos[0]} -v y=${pos[1]} -v z=${pos[2]} 'NR>9{print $x" "$y" "$z}' "$tmp_dump" > "$tmp_atom_coords"

        # Save the orientations. Make sure that the columns reserved for quaternions are not used to store velocities.
        if [[ ${#quat[@]} == 4 ]]; then 
           awk -v qw=${quat[0]} -v qx=${quat[1]} -v qy=${quat[2]} -v qz=${quat[3]} 'NR>9{print $qw" "$qx" "$qy" "$qz}' "$tmp_dump" > "$tmp_ellips_quat"
        fi

        # Save the velocities and, if present, angular momenta too.
        if [[ ${#vel[@]} == 3 && ${#angmom[@]} == 0 ]]; then 
           awk -v vx=${vel[0]} -v vy=${vel[1]} -v vz=${vel[2]} 'NR>9{print $1 " "$vx" "$vy" "$vz}' "$tmp_dump" > "$data_velocities"

        elif [[ ${#vel[@]} == 3 && ${#angmom[@]} == 3 ]]; then 
           awk -v vx=${vel[0]} -v vy=${vel[1]} -v vz=${vel[2]} -v ax=${angmom[0]} -v ay=${angmom[1]} -v az=${angmom[2]} 'NR>9{print $1 " "$vx" "$vy" "$vz" "$ax" "$ay" "$az}' "$tmp_dump" > "$data_velocities"
        fi
    
    elif [ "$A" = "-atomstyle" ] || [ "$A" = "-atom-style" ] || [ "$A" = "-atom_style" ]; then
        if [ "$i" -eq "$ARGC" ]; then
            echo "$SYNTAX_MSG" >&2
            exit 7
        fi
        i=$((i+1))
        eval A=\${ARGV${i}}
        if [ -z "$A" ]; then
            echo "$SYNTAX_MSG" >&2
            echo "-----------------------" >&2
            echo "" >&2
            echo "Error: The \"-atomstyle\" argument should be followed by an atom style." >&2
            echo "       (See the \"atom_style\" command in the LAMMPS documentation." >&2
            echo "        Note: hybrid atom styles are allowed but should be enclosed in quotes.)" >&2
            exit 8
        fi

        #echo "  (atom_style=\"$A\")" >&2
        #echo "" >&2

        ATOM_STYLE="$A"
        ATOM_STYLE_ARG="-atomstyle \"$A\""

        # (minor detail: $ATOM_STYLE_ARG should also be appended to TTREE_ARGS)
        if [ -z "$TTREE_ARGS" ]; then
            TTREE_ARGS="$ATOM_STYLE_ARG"
        else
            TTREE_ARGS="${TTREE_ARGS} $ATOM_STYLE_ARG"
        fi
	
    #else:  If the arguments are not understood in this script, then
    #       pass them on to "lttree.py"
    else
        A_FIRSTCHAR="$(echo $A| cut -c 1)"

        if [ "$A_FIRSTCHAR" = "\$" ]; then
            A="\\$A" # put an extra slash in front to prevent expansion later
        fi

        if [ -z "$TTREE_ARGS" ]; then
            TTREE_ARGS="\"$A\""
        else
            TTREE_ARGS="${TTREE_ARGS} \"$A\""
        fi
        # Check to see if this string ($A) ends in .lt or .LT
        # If so, then set the base name of the output files
        # to equal the base name of the .LT file being read.
        # (Being careful here.
        #  Sometimes the last argument is not the .lt or .LT file.
        #  Sometimes that file appears earlier in the argument list.
        #  I want to supply a default value.)
        #
        #   Note, in bash you can use:
        # if [ "${LAST_ARG/%.lt/}" -neq "$LAST_ARG" ]; then
        #     OUT_FILE_BASE="${LAST_ARG/%.lt/}"
        # But in the original bourn shell (sh), this does not work.
        # Instead we use a hack involving basename and dirname:

        if [ "$A_FIRSTCHAR" != "-" ]; then
            DN=`dirname "$A"`
            if [ "$DN" = "." ]; then
                DN=""
            else
                DN="${DN}/"
            fi

            BN=`basename "$A" .lt`
            if [ "${DN}${BN}" != "$A" ]; then
                OUT_FILE_BASE="$BN"
            else
                BN=`basename "$A" .LT`
                if [ "${DN}${BN}" != "$A" ]; then
                    OUT_FILE_BASE="$BN"
                fi
            fi
        fi
    fi
done


if [ -z "$ATOM_STYLE" ]; then
  #echo '########################################################' >&2
  #echo '##            WARNING: atom_style unspecified         ##' >&2
  #echo '##              Assuming atom_style = \"full\"          ##' >&2
  #echo '########################################################' >&2
  ATOM_STYLE="full"
fi



OUT_FILE_INPUT_SCRIPT="${OUT_FILE_BASE}.in"
OUT_FILE_INIT="${OUT_FILE_BASE}.in.init"
OUT_FILE_SETTINGS="${OUT_FILE_BASE}.in.settings"
OUT_FILE_DATA="${OUT_FILE_BASE}.data"
OUT_FILE_COORDS="${OUT_FILE_BASE}.in.coords"


if [ "$OUT_FILE_EXAMPLE_SCRIPT" == "" ]; then
    OUT_FILE_EXAMPLE_SCRIPT="run.in.EXAMPLE"  # default LAMMPS input script example
    if [ "$OUT_FILE_BASE" != "system" ]; then
        # For users who choose custom .LT file names,
        # the execution script should have the same custom name:
        OUT_FILE_EXAMPLE_SCRIPT=$OUT_FILE_INPUT_SCRIPT
    fi
fi

#  Now delete any old files lying around which were
#  created by previous invocations of moltemplate.
#
# In earlier versions of moltemplate, I was careful only to delete the most
# most common moltemplate-created files before running moltemplate again.
#rm -f "$OUT_FILE_EXAMPLE_SCRIPT" "$OUT_FILE_INIT" "$OUT_FILE_SETTINGS" "$OUT_FILE_DATA" "$OUT_FILE_COORDS"
# However, now I now delete all files beginning with "${OUT_FILE_BASE}.in.*"
# (except for those ending in .lt).
for f in "${OUT_FILE_BASE}.data" "${OUT_FILE_BASE}".in.*; do
  if [[ ! "$f" =~ .*\.lt ]]; then
    # Although it's unlikely somebody would name their moltemplate source file
    # "system.in.lt", just in case they do, we check to see if the file
    # name ends in ".lt" before deleting it.  Note: This also works:
    #if [ "$f" == `echo "$f" | sed  's|lt$||g'` ]; then ...
    # Delete the file:
    rm -f "$f"
  fi
done





# --------------------------------------------------------------------
# --- Now run ttree/lttree to generate the file fragments we need. ---
# ------  (Afterwards, we will paste these fragments together.)  -----
# --------------------------------------------------------------------




# If checking is not disabled, then first check for common spelling errors.

if [ -n "$LTTREE_CHECK_COMMAND" ]; then
    if ! eval $LTTREE_CHECK_COMMAND $TTREE_ARGS $LTTREE_CHECK_ARGS; then
        exit 1
    fi
fi

#   --- Run ttree. ---
#
# 3, 2, 1, ...

if ! eval $LTTREE_COMMAND $TTREE_ARGS; then
    exit 2
fi

# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# > Traceback (most recent call last):
# >   File "./lttree.py", line...
#            .-.
#           (0.0)
#         '=.|m|.='
#         .='`"``=.
#
# Hopefully this does not happen.
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

echo "" >&2



# Invoking the $LTTREE_COMMAND should generate a file ("ttree_assignments.txt")
# containing all of the "counter" variables (counters corresponding to
# atoms, bonds, atom-types, bond-types, etc...).
# Later, it will be convenient to create a version of this file
# which only contains static counter variables
# (ie. counter variables beginning with @, not $).
awk '/@/ {print $0}' < ttree_assignments.txt > ttree_assignments_static.txt



# Now count the number of atom-types, bond-types, angle-types, etc...

# Find the number of atom types, bond types, angle types, dihedra types, etc...
# The next line will count the number of times "@/atom" appears in the file:
NATOMTYPES=`awk 'BEGIN{n=0} /^@\/atom:/{n++}END{print n}' < ttree_assignments.txt`
NBONDTYPES=`awk 'BEGIN{n=0} /^@\/bond:/{n++}END{print n}' < ttree_assignments.txt`
NANGLETYPES=`awk 'BEGIN{n=0} /^@\/angle:/{n++}END{print n}' < ttree_assignments.txt`
NDIHEDRALTYPES=`awk 'BEGIN{n=0} /^@\/dihedral:/{n++}END{print n}' < ttree_assignments.txt`
NIMPROPERTYPES=`awk 'BEGIN{n=0} /^@\/improper:/{n++}END{print n}' < ttree_assignments.txt`

# The next line will count the number of times "@{/atom" appears in the file:
NATOMTYPES_SP=`awk 'BEGIN{n=0} /^@{\/atom:/{n++}END{print n}' < ttree_assignments.txt`
NBONDTYPES_SP=`awk 'BEGIN{n=0} /^@{\/bond:/{n++}END{print n}' < ttree_assignments.txt`
NANGLETYPES_SP=`awk 'BEGIN{n=0} /^@{\/angle:/{n++}END{print n}' < ttree_assignments.txt`
NDIHEDRALTYPES_SP=`awk 'BEGIN{n=0} /^@{\/dihedral:/{n++}END{print n}' < ttree_assignments.txt`
NIMPROPERTYPES_SP=`awk 'BEGIN{n=0} /^@{\/improper:/{n++}END{print n}' < ttree_assignments.txt`
# Add these numbers together to get the total number of times either
# "@{atom" or "@{/atom" appears in the ttree_assignments.txt file.
# This is the number of different atom types in the system.
NATOMTYPES=$((NATOMTYPES + NATOMTYPES_SP))
NBONDTYPES=$((NBONDTYPES + NBONDTYPES_SP))
NANGLETYPES=$((NANGLETYPES + NANGLETYPES_SP))
NDIHEDRALTYPES=$((NDIHEDRALTYPES + NDIHEDRALTYPES_SP))
NIMPROPERTYPES=$((NIMPROPERTYPES + NIMPROPERTYPES_SP))


# OMR: Moltemplate can handle the creation of valid input decks for LAMMPS even
# if the input LT file doesn't contain any atom. This could happen if the simulation
# has a rerun command, or a box is created in the simulation script.
#
#if [[ $NATOMTYPES -eq 0 ]]; then
    # Moltemplate can be used as a simple hierarchical template renderer
    # that knows nothing about LAMMPS.  In that case NATOMTYPES is undefined.
    # In that case, terminate now and do not try to interpret the files.
    #exit 0
#fi


# Attempt to remove any the DOS (windows) return-cairrage characters
# from inside the standard LAMMPS files generated by the user:
# (Users are free to put whatever weird characters they want in other
#  (custom) auxilliary files.  But not in the standard LAMMPS files.)


IFS=$CR
for file in $MOLTEMPLATE_TEMP_FILES; do
    if [ -e "$file" ]; then
        #dos2unix < "$file" > "$file.dos2unix"
        tr -d '\r' < "$file" > "$file.dos2unix"
        rm -f "$file" >/dev/null 2>&1 || true
        mv -f "$file.dos2unix" "$file"
    fi
done
IFS=$OIFS


if [ -s "${data_atoms}" ]; then
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicate_atoms.py" \
                                   < "${data_atoms}" \
                                   > "${data_atoms}.tmp"; then
        ERR_INTERNAL
    fi
    mv -f "${data_atoms}.tmp" "${data_atoms}"
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicate_atoms.py" \
                                   < "${data_atoms}.template" \
                                   > "${data_atoms}.tmp"; then
        ERR_INTERNAL
    fi
    mv -f "${data_atoms}.tmp" "${data_atoms}.template"

    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/renumber_DATA_first_column.py" \
                                       < "${data_atoms}" \
                                       > "${data_atoms}.tmp"; then
        ERR_INTERNAL
    fi
    mv -f "${data_atoms}.tmp" "${data_atoms}"
else
    if [[ $NATOMTYPES -eq 0 ]] && [[ -n "$LTTREE_CHECK_COMMAND" ]]; then
        echo "Error: There are no atoms in your system. Suggestions:" >&2
        echo "" >&2
        echo "       Your files should contain at least one" >&2
        echo "           write(\"${data_atoms}\")" >&2
        echo "       command.  These commands are typically located somewhere in" >&2
        echo "       one of the molecule object(s) you have defined." >&2
        echo "" >&2
        echo "       Make sure that you have the correct number of curly parenthesis {}." >&2
        echo "       (Extra \"}\" parenthesis can cause this error.)" >&2
        echo "" >&2
        echo "       This error also occurs if your input files lack \"new\" commands." >&2
        echo "       Once you have defined a type of molecule, you must create a copy" >&2
        echo "       of it using \"new\", if you want it to appear in your simulation." >&2
        echo "       See the moltemplate manual or online tutorials for examples." >&2
        echo "" >&2
        echo "       (This error also occurs if you instantiated an object using \"new\"" >&2
        echo "       which you thought was a molecule, but it is actually only a" >&2
        echo "       namespace, a force-field name or category containing only the" >&2
        echo "       definitions of other molecules, lacking any atoms of its own.)" >&2
        echo "" >&2
	echo "       To override this error, run Moltemplate using the \"-nocheck\" argument." >&2
        echo "" >&2
	exit 200
    fi
fi



# ---------------- Interactions By Type -----------------
# At the time of writing, bonded-interactions-by-atom-type were not
# understood by LAMMPS.  These features require auxilliary python scripts.
# These data sections must be processed before everything else (because
# they effect the other data sections, and the ttree_assignments.txt file.)
# -------------------------------------------------------

if [ -s "${data_bond_list}.template" ]; then

    if [ ! -s "$data_bonds_by_type" ]; then
        echo "Error: You have a \"Data Bond List\", section somewhere"
        echo "       without a \"Data Bonds By Type\" section to support it."
        echo "       (Did you mean to use \"Data Bonds\" instead?)"
        echo "Details:"
        echo "       Unlike the \"Data Bonds\" section, the \"Data Bond List\" section"
        echo "       allows the user to omit the bond types.  Instead moltemplate attempts"
        echo "       to infer the type of bond by considering the pair of atom types."
        echo "       However you must define a \"Data Bonds By Type\" section"
        echo "       to make this feature work (or use \"Data Bonds\" instead)."

        exit 15
    fi
    echo "Looking up bond types according to atom type" >&2
    #-- Generate a file containing bondid bondtype atomid1 atomid2 --
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/bonds_by_type.py" \
            -atom-style "$ATOM_STYLE" \
            -atoms "${data_atoms}.template" \
            -bond-list "${data_bond_list}.template" \
            -bondsbytype "${data_bonds_by_type}.template" \
            -prefix '$/bond:bytype' > gen_bonds.template.tmp; then
        exit 4
    fi

    # ---- cleanup: ----
    # ---- Create or re-build the "${data_bonds}.template" file ----
    # Instert these lines into the "${data_bonds}.template" file which includes
    # the newly generated interactions. (Note: these are in .template format)

    cp gen_bonds.template.tmp new_bonds.template.tmp
    if [ -s "${data_bonds}.template" ]; then
        # Then append existing "Bonds" to the end of the generated interactions
        # (Hopefully this way they will override those interactions.)
        cat "${data_bonds}.template" >> new_bonds.template.tmp
    fi
    mv -f new_bonds.template.tmp "${data_bonds}.template"



    #  ------  THE NEXT STEP IS NOT CURRENTLY NEEDED ------
    # All of the $bond variables have already been created, they just lack types
    # However we will need to do this if the user wants to omits the bond-ids.
    # In case I plan to allow the user to omit bond-ids, I leave this code here.
    #
    #echo "(Repairing ttree_assignments.txt file after bonds added.)" >&2
    #
    ## ---- Repair the ttree_assignments.txt file ----
    ## The next 2 lines extract the variable names from data_new.template.tmp
    ## and instert them into the appropriate place in ttree_assignments.txt
    ## (renumbering the relevant variable-assignments to avoid clashes).
    #if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_fix_ttree_assignments.py" \
    #      '/bond' gen_bonds.template.tmp \
    #      < ttree_assignments.txt \
    #      > ttree_assignments.tmp; then
    #    exit 5
    #fi
    #
    #echo "(Rendering ttree_assignments.tmp file after bonds added.)" >&2
    #mv -f ttree_assignments.tmp ttree_assignments.txt
    #  ----------------------------------------------------



    # ---- Re-build (render) the "$data_bonds" file ----
    # Now substitute these variable values (assignments) into the variable
    # names present in the .template file.  (We want to convert the file from
    # a .template format into an ordinary (numeric) LAMMPS data-section format.)

    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/ttree_render.py" \
           ttree_assignments.txt \
           < "${data_bonds}.template" \
           > "$data_bonds"; then
        exit 6
    fi

    echo "" >&2

    rm -f gen_bonds.template.tmp new_bonds.template.tmp

    echo "" >&2

fi






FILE_angles_by_type1=""
FILE_angles_by_type2=""
#for FILE in "$data_angles_by_type"*.template; do

IFS=$(echo -en "\n\b")
for FILE in `ls -v "$data_angles_by_type"*.template 2> /dev/null`; do

    if [ ! -s "$FILE" ] || [ ! -s "$data_bonds" ]; then
        break;  # This handles with the special cases that occur when
                # 1) There are no bonds in your system
                # 2) "$data_angles_by_type"*.template matches nothing
    fi

    echo "Generating 3-body angle interactions by atom/bond type" >&2

    # Extract the text between parenthesis (if present, empty-str otherwise)
    # Example: FILE="Data Angles By Type (gaff_angle.py)"
    SUBGRAPH_SCRIPT=`echo "$FILE" | awk '/\(.*\)/ {print $0}' | cut -d'(' -f2-| cut -d')' -f 1`
    # Example: (continued) SUBGRAPH_SCRIPT should equal "gaff_angle.py"

    # The user can also override this choice:
    if [ -n "$SUBGRAPH_SCRIPT_ANGLES" ]; then
        SUBGRAPH_SCRIPT="$SUBGRAPH_SCRIPT_ANGLES"
    elif [ -n "$SUBGRAPH_SCRIPT" ]; then
        SUBGRAPH_SCRIPT_ANGLES="$SUBGRAPH_SCRIPT"
    fi

    if [ -z "$SUBGRAPH_SCRIPT" ]; then
        SUBGRAPH_SCRIPT="nbody_Angles.py"
    else
        echo "(using the rules in \"$SUBGRAPH_SCRIPT\")" >&2
        # if [ ! -s "${PY_SCR_DIR}/nbody_alt_symmetry/$SUBGRAPH_SCRIPT" ]; then
        #     echo "Error: File \"$SUBGRAPH_SCRIPT\" not found." >&2
            #     echo "       It should be located in this directory:" >&2
        #     echo "       ${PY_SCR_DIR}/nbody_alt_symmetry/" >&2
        #     exit 4
        # fi
    fi

    FILE_angles_by_type2="$FILE_angles_by_type1"
    FILE_angles_by_type1="$FILE"

    #-- Generate a file containing the list of interactions on separate lines --
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_by_type.py" \
            -subgraph "${SUBGRAPH_SCRIPT}" \
            -section "Angles" \
            -sectionbytype "Angles By Type" \
            -atom-style "$ATOM_STYLE" \
            -atoms "${data_atoms}.template" \
            -bonds "${data_bonds}.template" \
            -nbodybytype "${FILE}" \
            $CHECKFF \
            -prefix '$/angle:bytype' > gen_angles.template.tmp; then
        exit 4
    fi

    # ---- cleanup: ----
    # ---- Re-build the "${data_angles}.template" file ----
    # Instert these lines into the "${data_angles}.template" file which includes
    # the newly generated interactions. (Note: these are in .template format)

    cp gen_angles.template.tmp new_angles.template.tmp
    if [ -s "${data_angles}.template" ]; then
        # Then append existing "Angles" to the end of the generated interactions
        # (Hopefully this way they will override those interactions.)
        cat "${data_angles}.template" >> new_angles.template.tmp
    fi
    mv -f new_angles.template.tmp "${data_angles}.template"

    echo "(Repairing ttree_assignments.txt file after angles added.)" >&2

    # ---- Repair the ttree_assignments.txt file ----
    # The next 2 lines extract the variable names from data_new.template.tmp
    # and instert them into the appropriate place in ttree_assignments.txt
    # (renumbering the relevant variable-assignments to avoid clashes).
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_fix_ttree_assignments.py" \
          '/angle' gen_angles.template.tmp \
          < ttree_assignments.txt \
          > ttree_assignments.tmp; then
        exit 5
    fi

    echo "(Rendering ttree_assignments.tmp file after angles added.)" >&2

    # ---- Re-build (render) the "$data_angles" file ----
    # Now substitute these variable values (assignments) into the variable
    # names present in the .template file.  (We want to convert the file from
    # a .template format into an ordinary (numeric) LAMMPS data-section format.)
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/ttree_render.py" \
           ttree_assignments.tmp \
           < "${data_angles}.template" \
           > "$data_angles"; then
        exit 6
    fi
    echo "" >&2

    mv -f ttree_assignments.tmp ttree_assignments.txt
    rm -f gen_angles.template.tmp new_angles.template.tmp
done
IFS=$OIFS



FILE_dihedrals_by_type1=""
FILE_dihedrals_by_type2=""
#for FILE in "$data_dihedrals_by_type"*.template; do
IFS=$(echo -en "\n\b")
for FILE in `ls -v "$data_dihedrals_by_type"*.template 2> /dev/null`; do

    if [ ! -s "$FILE" ] || [ ! -s "$data_bonds" ]; then
        break;  # This handles with the special cases that occur when
                # 1) There are no bonds in your system
                # 2) "$data_dihedrals_by_type"*.template matches nothing
    fi

    echo "Generating 4-body dihedral interactions by atom/bond type" >&2

    # Extract the text between parenthesis (if present, empty-str otherwise)
    # Example: FILE="Data Dihedrals By Type (gaff_dih.py)"
    SUBGRAPH_SCRIPT=`echo "$FILE" | awk '/\(.*\)/ {print $0}' | cut -d'(' -f2-| cut -d')' -f 1`
    # Example: (continued) SUBGRAPH_SCRIPT should equal "gaff_dih.py"

    # The user can also override this choice:
    if [ -n "$SUBGRAPH_SCRIPT_DIHEDRALS" ]; then
        SUBGRAPH_SCRIPT="$SUBGRAPH_SCRIPT_DIHEDRALS"
    elif [ -n "$SUBGRAPH_SCRIPT" ]; then
        SUBGRAPH_SCRIPT_DIHEDRALS="$SUBGRAPH_SCRIPT"
    fi

    if [ -z "$SUBGRAPH_SCRIPT" ]; then
        SUBGRAPH_SCRIPT="nbody_Dihedrals.py"
    else
        echo "(using the rules in \"$SUBGRAPH_SCRIPT\")" >&2
        # if [ ! -s "${PY_SCR_DIR}/nbody_alt_symmetry/$SUBGRAPH_SCRIPT" ]; then
        #     echo "Error: File \"$SUBGRAPH_SCRIPT\" not found." >&2
            # echo "       It should be located in this directory:" >&2
        #     echo "       ${PY_SCR_DIR}/nbody_alt_symmetry/" >&2
        #     exit 4
        # fi
    fi

    FILE_dihedrals_by_type2="$FILE_dihedrals_by_type1"
    FILE_dihedrals_by_type1="$FILE"

    #-- Generate a file containing the list of interactions on separate lines --
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_by_type.py" \
            -subgraph "${SUBGRAPH_SCRIPT}" \
            -section "Dihedrals" \
            -sectionbytype "Dihedrals By Type" \
            -atom-style "$ATOM_STYLE" \
            -atoms "${data_atoms}.template" \
            -bonds "${data_bonds}.template" \
            -nbodybytype "${FILE}" \
            $CHECKFF \
            -prefix '$/dihedral:bytype' > gen_dihedrals.template.tmp; then
        exit 4
    fi

    # ---- cleanup: ----
    # ---- Re-build the "${data_dihedrals}.template" file ----
    # Instert these lines into the "${data_dihedrals}.template" file which includes
    # the newly generated interactions. (Note: these are in .template format)

    cp gen_dihedrals.template.tmp new_dihedrals.template.tmp
    if [ -s "${data_dihedrals}.template" ]; then
        # Then append existing "Dihedrals" to the end of the generated interactions
        # (Hopefully this way they will override those interactions.)
        cat "${data_dihedrals}.template" >> new_dihedrals.template.tmp
    fi
    mv -f new_dihedrals.template.tmp "${data_dihedrals}.template"

    echo "(Repairing ttree_assignments.txt file after dihedrals added.)" >&2

    # ---- Repair the ttree_assignments.txt file ----
    # The next 2 lines extract the variable names from data_new.template.tmp
    # and instert them into the appropriate place in ttree_assignments.txt
    # (renumbering the relevant variable-assignments to avoid clashes).
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_fix_ttree_assignments.py" \
          '/dihedral' gen_dihedrals.template.tmp \
          < ttree_assignments.txt \
          > ttree_assignments.tmp; then
        exit 5
    fi

    echo "(Rendering ttree_assignments.tmp file after dihedrals added.)" >&2

    # ---- Re-build (render) the "$data_dihedrals" file ----
    # Now substitute these variable values (assignments) into the variable
    # names present in the .template file.  (We want to convert the file from
    # a .template format into an ordinary (numeric) LAMMPS data-section format.)
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/ttree_render.py" \
           ttree_assignments.tmp \
           < "${data_dihedrals}.template" \
           > "$data_dihedrals"; then
        exit 6
    fi
    echo "" >&2

    mv -f ttree_assignments.tmp ttree_assignments.txt
    rm -f gen_dihedrals.template.tmp new_dihedrals.template.tmp
done
IFS=$OIFS





FILE_impropers_by_type1=""
FILE_impropers_by_type2=""
#for FILE in "$data_impropers_by_type"*.template; do
IFS=$(echo -en "\n\b")
for FILE in `ls -v "$data_impropers_by_type"*.template 2> /dev/null`; do

    if [ ! -s "$FILE" ] || [ ! -s "$data_bonds" ]; then
        break;  # This handles with the special cases that occur when
                # 1) There are no bonds in your system
                # 2) "$data_impropers_by_type"*.template matches nothing
    fi

    echo "Generating 4-body improper interactions by atom/bond type" >&2

    # Extract the text between parenthesis (if present, empty-str otherwise)
    # Example: FILE="Data Impropers By Type (gaff_impr.py)"
    SUBGRAPH_SCRIPT=`echo "$FILE" | awk '/\(.*\)/ {print $0}' | cut -d'(' -f2-| cut -d')' -f 1`
    # Example: (continued) SUBGRAPH_SCRIPT should equal "gaff_impr.py"

    # The user can also override this choice:
    if [ -n "$SUBGRAPH_SCRIPT_IMPROPERS" ]; then
        SUBGRAPH_SCRIPT="$SUBGRAPH_SCRIPT_IMPROPERS"
    elif [ -n "$SUBGRAPH_SCRIPT" ]; then
        SUBGRAPH_SCRIPT_IMPROPERS="$SUBGRAPH_SCRIPT"
    fi

    if [ -z "$SUBGRAPH_SCRIPT" ]; then
        SUBGRAPH_SCRIPT="nbody_Impropers.py"
    else
        echo "(using the rules in \"$SUBGRAPH_SCRIPT\")" >&2
        # if [ ! -s "${PY_SCR_DIR}/nbody_alt_symmetry/$SUBGRAPH_SCRIPT" ]; then
        #     echo "Error: File \"$SUBGRAPH_SCRIPT\" not found." >&2
            # echo "       It should be located in this directory:" >&2
        #     echo "       ${PY_SCR_DIR}/nbody_alt_symmetry/" >&2
        #     exit 4
        # fi
    fi

    FILE_impropers_by_type2="$FILE_impropers_by_type1"
    FILE_impropers_by_type1="$FILE"

    #-- Generate a file containing the list of interactions on separate lines --
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_by_type.py" \
            -subgraph "${SUBGRAPH_SCRIPT}" \
            -section "Impropers" \
            -sectionbytype "Impropers By Type" \
            -atom-style "$ATOM_STYLE" \
            -atoms "${data_atoms}.template" \
            -bonds "${data_bonds}.template" \
            -nbodybytype "${FILE}" \
            -prefix '$/improper:bytype' > gen_impropers.template.tmp; then
        exit 4
    fi

    # ---- cleanup: ----
    # ---- Re-build the "${data_impropers}.template" file ----
    # Instert these lines into the "${data_impropers}.template" file which includes
    # the newly generated interactions. (Note: these are in .template format)

    cp gen_impropers.template.tmp new_impropers.template.tmp
    if [ -s "${data_impropers}.template" ]; then
        # Then append existing "Impropers" to the end of the generated interactions
        # (Hopefully this way they will override those interactions.)
        cat "${data_impropers}.template" >> new_impropers.template.tmp
    fi
    mv -f new_impropers.template.tmp "${data_impropers}.template"

    echo "(Repairing ttree_assignments.txt file after impropers added.)" >&2

    # ---- Repair the ttree_assignments.txt file ----
    # The next 2 lines extract the variable names from data_new.template.tmp
    # and instert them into the appropriate place in ttree_assignments.txt
    # (renumbering the relevant variable-assignments to avoid clashes).
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_fix_ttree_assignments.py" \
          '/improper' gen_impropers.template.tmp \
          < ttree_assignments.txt \
          > ttree_assignments.tmp; then
        exit 5
    fi

    echo "(Rendering ttree_assignments.tmp file after impropers added.)" >&2

    # ---- Re-build (render) the "$data_impropers" file ----
    # Now substitute these variable values (assignments) into the variable
    # names present in the .template file.  (We want to convert the file from
    # a .template format into an ordinary (numeric) LAMMPS data-section format.)
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/ttree_render.py" \
           ttree_assignments.tmp \
           < "${data_impropers}.template" \
           > "$data_impropers"; then
        exit 6
    fi
    echo "" >&2

    mv -f ttree_assignments.tmp ttree_assignments.txt
    rm -f gen_impropers.template.tmp new_impropers.template.tmp
done
IFS=$OIFS



# Find all the files created by lttree.py containing lines beginning with 
# "pair_coeff", "bond_coeff", "angle_coeff", "dihedral_coeff", "improper_coeff".
# These files need to be processed further.


OUT_FILES_WITH_COEFF_COMMANDS=""
for f in *.template; do
    HAS_COEFF_COMMANDS=`awk 'BEGIN{has_coeffs=0} {if ((NF>0)&&(($1=="pair_coeff")||($1=="bond_coeff")||($1=="angle_coeff")||($1=="dihedral_coeff")||($1=="improper_coeff"))) has_coeffs=1} END{print has_coeffs}' < "$f"`
    if [ "$HAS_COEFF_COMMANDS" -eq "1" ]; then
        if [ -n "$OUT_FILES_WITH_COEFF_COMMANDS" ]; then
            printf -v OUT_FILES_WITH_COEFF_COMMANDS "${OUT_FILES_WITH_COEFF_COMMANDS}\n$f"
        else
            OUT_FILES_WITH_COEFF_COMMANDS="$f"
        fi
    fi
done

# Deal with wildcard characters ('*', '?') in "_coeff" commands
# appearing in any LAMMPS input scripts generated by moltemplate.
# Replace them with explicit variable names.  Do this before rendering
#if [ -s "${in_settings}.template" ]; then
echo "expanding wildcards in \"_coeff\" commands" >&2

IFS=$CR
for file_name in $OUT_FILES_WITH_COEFF_COMMANDS; do
    #echo "  searching for \"_coeff\" commands in file \"$file_name\"" >&2
    # file contains both _coeff commands and wildcards *,? on the same line:
    if ! awk '{if (match($1,/'_coeff/') && match($0,/'[*,?]/')) exit 1}' < "$file_name"; then

        echo "  expanding wildcards in \"_coeff\" commands in \"$file_name\"">&2
        if ! eval $PYTHON_COMMAND "${PY_SCR_DIR}/postprocess_coeffs.py" ttree_assignments_static.txt < "$file_name" > "${file_name}.tmp"; then
            ERR_INTERNAL
        fi

        mv -f "${file_name}.tmp" "$file_name"
        # Now reassign integers to these variables
        bn=`basename "$file_name" .template`
        if grep '\$' < "$file_name" > /dev/null; then
            # If "$file_name" contains the $ character, then there might be
            # $-style counter variables somewhere in the file.  In that case
            # we must consider all the variables in the entire
            # "ttree_assignments.txt" file, when substituting numbers for
            # variables.  This is a very big file and can take a while to read
            # so we don't do it unless it's necessary.
            if ! $PYTHON_COMMAND "${PY_SCR_DIR}/ttree_render.py" \
                 ttree_assignments.txt \
                 < "$file_name" \
                 > "$bn"; then
                exit 6
            fi
        else
            # If "$file_name" does not contains the $ character, then we can
            # ignore all of the $-style variables in the "ttree_assignments.txt"
            # file, which will save a lot of time.  In that case we can use
            # "ttree_assignments_static.txt" instead which omits those lines
            # of text and is much faster to parse as a result.
            if ! $PYTHON_COMMAND "${PY_SCR_DIR}/ttree_render.py" \
                 ttree_assignments_static.txt \
                 < "$file_name" \
                 > "$bn"; then
                exit 6
            fi
        fi
    fi
done
#fi
IFS=$OIFS

if [ -n "$LTTREE_POSTPROCESS_COMMAND" ]; then
    echo "" >&2
    if ! eval $LTTREE_POSTPROCESS_COMMAND $TTREE_ARGS; then
        exit 3
    fi
    echo "" >&2
fi


# -------------------------------------------------------
# If present, then remove duplicate bonds, angles, dihedrals, and impropers
# (unless overridden by the user).
# -------------------------------------------------------


if [ -s "${data_masses}" ]; then
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicate_atoms.py" \
                                   < "${data_masses}" \
                                   > "${data_masses}.tmp"; then
        ERR_INTERNAL
    fi
    mv -f "${data_masses}.tmp" "${data_masses}"
fi


if [ -s "${data_bonds}" ]; then
    SUBGRAPH_SCRIPT="nbody_Bonds.py"
    if [ -n "$SUBGRAPH_SCRIPT_BONDS" ]; then
        SUBGRAPH_SCRIPT="$SUBGRAPH_SCRIPT_BONDS"
    fi
    if [ ! -z $REMOVE_DUPLICATE_BONDS ]; then
        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_reorder_atoms.py" \
                             Bonds \
                             "$SUBGRAPH_SCRIPT" \
                             < "${data_bonds}" \
                             > "${data_bonds}.tmp"; then
            ERR_INTERNAL
        fi
        cp -f "${data_bonds}.tmp" "${data_bonds}"
        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 2 \
                             < "${data_bonds}" \
                             > "${data_bonds}.tmp"; then
            ERR_INTERNAL
        fi
        mv "${data_bonds}.tmp" "${data_bonds}"
	if [ ! -z "$REPORT_DUPLICATES" ]; then
          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 2 Bond warning_duplicate_bonds.txt "$RD_ID_FILTER" "$RD_TYPE_FILTER" \
                               < "${data_bonds}.template" \
                               > "${data_bonds}.tmp"; then
              ERR_INTERNAL
          fi
	else
          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 2 \
                               < "${data_bonds}.template" \
                               > "${data_bonds}.tmp"; then
              ERR_INTERNAL
          fi
	fi
        mv "${data_bonds}.tmp" "${data_bonds}.template"
    fi
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/renumber_DATA_first_column.py" \
                         < "${data_bonds}" \
                         > "${data_bonds}.tmp"; then
        ERR_INTERNAL
    fi
    mv -f "${data_bonds}.tmp" "${data_bonds}"
fi




if [ -s "${data_angles}" ]; then
    SUBGRAPH_SCRIPT="nbody_Angles.py"
    if [ -n "$SUBGRAPH_SCRIPT_ANGLES" ]; then
        SUBGRAPH_SCRIPT="$SUBGRAPH_SCRIPT_ANGLES"
    fi
    if [ ! -z $REMOVE_DUPLICATE_ANGLES ]; then
        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_reorder_atoms.py" \
                             Angles \
                             "$SUBGRAPH_SCRIPT" \
                             < "${data_angles}" \
                             > "${data_angles}.tmp"; then
            ERR_INTERNAL
        fi
        cp -f "${data_angles}.tmp" "${data_angles}"
        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 3 \
                             < "${data_angles}" \
                             > "${data_angles}.tmp"; then
            ERR_INTERNAL
        fi
        mv "${data_angles}.tmp" "${data_angles}"
	if [ ! -z "$REPORT_DUPLICATES" ]; then
          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 3 Angle warning_duplicate_angles.txt "$RD_ID_FILTER" "$RD_TYPE_FILTER" \
                               < "${data_angles}.template" \
                               > "${data_angles}.tmp"; then
              ERR_INTERNAL
          fi
	else
	  if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 3 \
                               < "${data_angles}.template" \
                               > "${data_angles}.tmp"; then
              ERR_INTERNAL
          fi
        fi
        mv "${data_angles}.tmp" "${data_angles}".template
    fi
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/renumber_DATA_first_column.py" \
                         < "${data_angles}" \
                         > "${data_angles}.tmp"; then
        ERR_INTERNAL
    fi
    mv -f "${data_angles}.tmp" "${data_angles}"

    if  [ ! -z $FILE_angles_by_type2 ]; then
        MSG_MULTIPLE_ANGLE_RULES=$(cat <<EOF
#############################################################################
WARNING:
  It appears as though multiple conflicting rules were used to generate
ANGLE interactions.  (This can occur when combining molecules built with
different force-field rules).  In your case, you are using rules defined here:
   "$FILE_angles_by_type2"
   "$FILE_angles_by_type1"
   (Files ending in .py are located here:
    $PY_SCR_DIR/nbody_alt_symmetry/)
If the molecules built using these two different force-field settings are not
connected, AND if you do NOT override force-field angles with explicitly
defined angles, then you can probably ignore this warning message.  Otherwise
please check the list of angle interactions to make sure they are correct!
(It might help to build a much smaller system using the same molecule types.)
#############################################################################

EOF
)
        echo "$MSG_MULTIPLE_ANGLE_RULES" >&2
    fi
fi





if [ -s "${data_dihedrals}" ]; then
    SUBGRAPH_SCRIPT="nbody_Dihedrals.py"
    if [ -n "$SUBGRAPH_SCRIPT_DIHEDRALS" ]; then
        SUBGRAPH_SCRIPT="$SUBGRAPH_SCRIPT_DIHEDRALS"
    fi
    if [ ! -z $REMOVE_DUPLICATE_DIHEDRALS ]; then
        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_reorder_atoms.py" \
                             Dihedrals \
                             "$SUBGRAPH_SCRIPT" \
                             < "${data_dihedrals}" \
                             > "${data_dihedrals}.tmp"; then
            ERR_INTERNAL
        fi
        cp -f "${data_dihedrals}.tmp" "${data_dihedrals}"
        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 \
                             < "${data_dihedrals}" \
                             > "${data_dihedrals}.tmp"; then
            ERR_INTERNAL
        fi
        mv "${data_dihedrals}.tmp" "${data_dihedrals}"
	if [ ! -z "$REPORT_DUPLICATES" ]; then
          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 Dihedral warning_duplicate_dihedrals.txt "$RD_ID_FILTER" "$RD_TYPE_FILTER" \
                               < "${data_dihedrals}.template" \
                               > "${data_dihedrals}.tmp"; then
              ERR_INTERNAL
          fi
	else
          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 \
                               < "${data_dihedrals}.template" \
                               > "${data_dihedrals}.tmp"; then
              ERR_INTERNAL
          fi
	fi
        mv "${data_dihedrals}.tmp" "${data_dihedrals}.template"
    fi
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/renumber_DATA_first_column.py" \
                         < "${data_dihedrals}" \
                         > "${data_dihedrals}.tmp"; then
        ERR_INTERNAL
    fi
    mv -f "${data_dihedrals}.tmp" "${data_dihedrals}"

    if  [ ! -z $FILE_dihedrals_by_type2 ]; then
        MSG_MULTIPLE_DIHEDRAL_RULES=$(cat <<EOF
#############################################################################
WARNING:
  It appears as though multiple conflicting rules were used to generate
DIHEDRAL interactions.  (This can occur when combining molecules built with
different force-field rules).  In your case, you are using rules defined here:
   "$FILE_dihedrals_by_type2"
   "$FILE_dihedrals_by_type1"
   (Files ending in .py are located here:
    $PY_SCR_DIR/nbody_alt_symmetry/)
If the molecules built using these two different force-field settings are not
connected, AND if you do NOT override force-field dihedrals with explicitly
defined dihedrals, then you can probably ignore this warning message.  Otherwise
please check the list of dihedral interactions to make sure they are correct!
(It might help to build a much smaller system using the same molecule types.)
#############################################################################

EOF
)
        echo "$MSG_MULTIPLE_DIHEDRAL_RULES" >&2
    fi
fi





if [ -s "${data_impropers}" ]; then
    SUBGRAPH_SCRIPT="nbody_Impropers.py"
    if [ -n "$SUBGRAPH_SCRIPT_IMPROPERS" ]; then
        SUBGRAPH_SCRIPT="$SUBGRAPH_SCRIPT_IMPROPERS"
    fi
    if [ ! -z $REMOVE_DUPLICATE_IMPROPERS ]; then
        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/nbody_reorder_atoms.py" \
                             Impropers \
                             "$SUBGRAPH_SCRIPT" \
                             < "${data_impropers}" \
                             > "${data_impropers}.tmp"; then
            ERR_INTERNAL
        fi
        cp -f "${data_impropers}.tmp" "${data_impropers}"
        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 \
                             < "${data_impropers}" \
                             > "${data_impropers}.tmp"; then
            ERR_INTERNAL
        fi
        mv "${data_impropers}.tmp" "${data_impropers}"
	if [ ! -z "$REPORT_DUPLICATES" ]; then
          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 Improper warning_duplicate_impropers.txt "$RD_ID_FILTER" "$RD_TYPE_FILTER" \
                               < "${data_impropers}.template" \
                               > "${data_impropers}.tmp"; then
              ERR_INTERNAL
          fi
	else
          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 \
                               < "${data_impropers}.template" \
                               > "${data_impropers}.tmp"; then
              ERR_INTERNAL
          fi
	fi
        mv "${data_impropers}.tmp" "${data_impropers}.template"
    fi
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/renumber_DATA_first_column.py" \
                         < "${data_impropers}" \
                         > "${data_impropers}.tmp"; then
        ERR_INTERNAL
    fi
    mv -f "${data_impropers}.tmp" "${data_impropers}"

    if  [ ! -z $FILE_impropers_by_type2 ]; then
        MSG_MULTIPLE_IMPROPER_RULES=$(cat <<EOF
#############################################################################
WARNING:
  It appears as though multiple conflicting rules were used to generate
IMPROPER interactions.  (This can occur when combining molecules built with
different force-field rules.)  In your case, you are using rules defined here:
   "$FILE_impropers_by_type2"
   "$FILE_impropers_by_type1"
   (Files ending in .py are located here:
    $PY_SCR_DIR/nbody_alt_symmetry/)
If the molecules built using these two different force-field settings are not
connected, AND if you do NOT override force-field imrpopers with explicitly
defined impropers, then you can probably ignore this warning message.  Otherwise
please check the list of improper interactions to make sure they are correct!
(It might help to build a much smaller system using the same molecule types.)
#############################################################################

EOF
)
        echo "$MSG_MULTIPLE_IMPROPER_RULES" >&2
    fi

fi





# ------------------ Charge By Bond ----------------------
# Assign atom partial charges according to who they are bonded to

if [ -s "$data_charge_by_bond" ]; then
    echo "Looking up partial charge contributions from bonds" >&2
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/charge_by_bond.py" \
         -atom-style "$ATOM_STYLE" \
         -atoms "${data_atoms}.template" \
         -bonds "${data_bonds}.template" \
         -bond-list "${data_bond_list}.template" \
         -chargebybond "${data_charge_by_bond}.template" \
         > gen_charges.template.tmp; then
        exit 4
    fi

    # ---- cleanup: ----
    # ---- Create or re-build the "${in_charges}.template" file ----
    # Instert these lines into the "${in_charges}.template" file which includes
    # the newly generated interactions. (Note: these are in .template format)

    cp gen_charges.template.tmp new_charges.template.tmp
    if [ -s "${in_charges}.template" ]; then
        # Then append existing "Bonds" to the end of the generated interactions
        # (Hopefully this way they will override those interactions.)
        cat "${in_charges}.template" >> new_charges.template.tmp
    fi
    mv -f new_charges.template.tmp "${in_charges}.template"

    # ---- Re-build (render) the "$in_charges" file ----
    # Now substitute these variable values (assignments) into the variable
    # names present in the .template file.  (We want to convert the file from
    # a .template format into an ordinary (numeric) LAMMPS data-section format.)
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/ttree_render.py" \
         ttree_assignments.txt \
         < "${in_charges}.template" \
         >> "${in_charges}"; then
        exit 6
    fi
    echo "" >&2

    rm -f gen_charges.template.tmp new_charges.template.tmp

    echo "" >&2
fi



# -------------------------------------------------------

rm -f "$OUT_FILE_DATA"

#NATOMS=`awk 'BEGIN{n=0} /^\\\$\/atom:/{n++}END{print n}' < ttree_assignments.txt`
#NBONDS=`awk 'BEGIN{n=0} /^\\\$\/bond:/{n++}END{print n}' < ttree_assignments.txt`
#NANGLES=`awk 'BEGIN{n=0} /^\\\$\/angle:/{n++}END{print n}' < ttree_assignments.txt`
#NDIHEDRALS=`awk 'BEGIN{n=0} /^\\\$\/dihedral:/{n++}END{print n}' < ttree_assignments.txt`
#NIMPROPERS=`awk 'BEGIN{n=0} /^\\\$\/improper:/{n++}END{print n}' < ttree_assignments.txt`
#
#NATOMS_SP=`awk 'BEGIN{n=0} /^\\\${\/atom:/{n++}END{print n}' < ttree_assignments.txt`
#NBONDS_SP=`awk 'BEGIN{n=0} /^\\\${\/bond:/{n++}END{print n}' < ttree_assignments.txt`
#NANGLES_SP=`awk 'BEGIN{n=0} /^\\\${\/angle:/{n++}END{print n}' < ttree_assignments.txt`
#NDIHEDRALS_SP=`awk 'BEGIN{n=0} /^\\\${\/dihedral:/{n++}END{print n}' < ttree_assignments.txt`
#NIMPROPERS_SP=`awk 'BEGIN{n=0} /^\\\${\/improper:/{n++}END{print n}' < ttree_assignments.txt`
#
#NATOMS=$((NATOMS + NATOMS_SP))
#NBONDS=$((NBONDS + NBONDS_SP))
#NANGLES=$((NANGLES + NANGLES_SP))
#NDIHEDRALS=$((NDIHEDRALS + NDIHEDRALS_SP))
#NIMPROPERS=$((NIMPROPERS + NIMPROPERS_SP))

NATOMS="0"
NBONDS="0"
NANGLES="0"
NDIHEDRALS="0"
NIMPROPERS="0"
NCMAP="0"

if [ -s "${data_atoms}" ]; then
  NATOMS=`awk 'END{print NR}' < "${data_atoms}"`
fi
if [ -s "${data_ellipsoids}" ]; then
  NELLIPSOIDS=`awk 'END{print NR}' < "${data_ellipsoids}"`
fi
if [ -s "${data_bonds}" ]; then
  NBONDS=`awk 'END{print NR}' < "${data_bonds}"`
fi
if [ -s "${data_angles}" ]; then
  NANGLES=`awk 'END{print NR}' < "${data_angles}"`
fi
if [ -s "${data_dihedrals}" ]; then
  NDIHEDRALS=`awk 'END{print NR}' < "${data_dihedrals}"`
fi
if [ -s "${data_impropers}" ]; then
  NIMPROPERS=`awk 'END{print NR}' < "${data_impropers}"`
fi
if [ -s "${data_cmap}" ]; then
  NCMAP=`awk 'END{print NR}' < "${data_cmap}"`
fi


echo "LAMMPS Description" > "$OUT_FILE_DATA"
echo "" >> "$OUT_FILE_DATA"
echo "     $NATOMS  atoms" >> "$OUT_FILE_DATA"
if [ -n "$NELLIPSOIDS" ]; then
    echo "     $NELLIPSOIDS  ellipsoids" >> "$OUT_FILE_DATA"
fi
if [ -n "$NBONDS" ] || [ -n "$NBONDTYPES" ]; then
    echo "     $NBONDS  bonds" >> "$OUT_FILE_DATA"
fi
if [ -n "$NANGLES" ] || [ -n "$NANGLETYPES" ]; then
    echo "     $NANGLES  angles" >> "$OUT_FILE_DATA"
fi
if [ -n "$NDIHEDRALS" ] || [ -n "$NDIHEDRALTYPES" ]; then
    echo "     $NDIHEDRALS  dihedrals" >> "$OUT_FILE_DATA"
fi
if [ -n "$NIMPROPERS" ] || [ -n "$NIMPROPERTYPES" ]; then
    echo "     $NIMPROPERS  impropers" >> "$OUT_FILE_DATA"
fi
if [ -n "$NCMAP" ] && [ "$NCMAP" -gt "0" ]; then
   echo "     $NCMAP  crossterms" >> "$OUT_FILE_DATA"
fi

echo "" >> "$OUT_FILE_DATA"
echo "     $NATOMTYPES  atom types" >> "$OUT_FILE_DATA"

if [ -n "$NBONDTYPES" ]; then
    echo "     $NBONDTYPES  bond types" >> "$OUT_FILE_DATA"
fi
if [ -n "$NANGLETYPES" ]; then
    echo "     $NANGLETYPES  angle types" >> "$OUT_FILE_DATA"
fi
if [ -n "$NDIHEDRALTYPES" ]; then
    echo "     $NDIHEDRALTYPES  dihedral types" >> "$OUT_FILE_DATA"
fi
if [ -n "$NIMPROPERTYPES" ]; then
    echo "     $NIMPROPERTYPES  improper types" >> "$OUT_FILE_DATA"
fi
echo "" >> "$OUT_FILE_DATA"




if [ -s "$data_header" ]; then
  cat "$data_header" >> "$OUT_FILE_DATA"
  echo "" >> "$OUT_FILE_DATA"
fi




# --- PERIODIC BOUNDARY CONDITIONS ---

# Note: If there is a "$data_boundary" file present, it overrides any settings
#       which may have been stored in a pdb file or other external file.

if [ -s "$data_pbc" ] && [ ! -s "$data_boundary" ]; then
    mv -f "$data_pbc" "$data_boundary"
    echo "WARNING: write_once(\"$data_pbc\") is depreciated" >&2
    echo "     Use write_once(\"$data_boundary\") instead" >&2
fi

if [ -s "$data_boundary" ]; then
    # Copy the boundary conditions from the "$data_boundary" file.
    # Don't assume there is only one line containing "xlo xhi", for example.
    # It's possible the user wrote the boundary conditions multiple times.
    # As always, the most recent setting overrides the earlier settings.
    BOXSIZE_MINX=`tr -d '\015' < "$data_boundary" | awk '{if ($3=="xlo") {xlo=$1}} END{print xlo}'`
    BOXSIZE_MAXX=`tr -d '\015' < "$data_boundary" | awk '{if ($4=="xhi") {xhi=$2}} END{print xhi}'`

    BOXSIZE_MINY=`tr -d '\015' < "$data_boundary" | awk '{if ($3=="ylo") {ylo=$1}} END{print ylo}'`
    BOXSIZE_MAXY=`tr -d '\015' < "$data_boundary" | awk '{if ($4=="yhi") {yhi=$2}} END{print yhi}'`

    BOXSIZE_MINZ=`tr -d '\015' < "$data_boundary" | awk '{if ($3=="zlo") {zlo=$1}} END{print zlo}'`
    BOXSIZE_MAXZ=`tr -d '\015' < "$data_boundary" | awk '{if ($4=="zhi") {zhi=$2}} END{print zhi}'`

    if [ -z "$BOXSIZE_MINX" ] || [ -z "$BOXSIZE_MAXX" ]; then
        echo "Error: Problem with box boundary format (\"xlo xhi\") in \"$data_boundary\"" >&2
        exit 12
    fi
    if [ -z "$BOXSIZE_MINY" ] || [ -z "$BOXSIZE_MAXY" ]; then
        echo "Error: Problem with box boundary format (\"ylo yhi\") in \"$data_boundary\"" >&2
        exit 12
    fi
    if [ -z "$BOXSIZE_MINZ" ] || [ -z "$BOXSIZE_MAXZ" ]; then
        echo "Error: Problem with box boundary format (\"zlo zhi\") in \"$data_boundary\"" >&2
        exit 12
    fi

    BOXSIZE_XY=`awk '{if ($4=="xy") {xy=$1}} END{print xy}' < "$data_boundary"`
    BOXSIZE_XZ=`awk '{if ($5=="xz") {xz=$2}} END{print xz}' < "$data_boundary"`
    BOXSIZE_YZ=`awk '{if ($6=="yz") {yz=$3}} END{print yz}' < "$data_boundary"`

    if [ -n "$BOXSIZE_XY" ] || [ -n "$BOXSIZE_XZ" ] || [ -n "$BOXSIZE_YZ" ]; then
        if [ -n "$BOXSIZE_XY" ] && [ -n "$BOXSIZE_XZ" ] && [ -n "$BOXSIZE_YZ" ]; then
            #echo "triclinic_parameters: XY XZ YZ = $BOXSIZE_XY $BOXSIZE_XZ $BOXSIZE_YZ" >&2
            TRICLINIC="true"
        else
            echo "Error: Problem with triclinic format (\"xy xz yz\") in \"$data_boundary\"" >&2
            exit 13
        fi
    fi
fi




if [ -z "$BOXSIZE_MINX" ] || [ -z "$BOXSIZE_MAXX" ] || [ -z "$BOXSIZE_MINY" ] || [ -z "$BOXSIZE_MAXY" ] || [ -z "$BOXSIZE_MINZ" ] || [ -z "$BOXSIZE_MAXZ" ]; then
    echo "Periodic boundary conditions unspecified. Attempting to generate automatically." >&2

    # By default, disable triclinic
    BOXSIZE_XY=""
    BOXSIZE_XZ=""
    BOXSIZE_YZ=""
    TRICLINIC=""

    if [ -s "$tmp_atom_coords" ]; then
        # Estimate the minimimum, maximum x,y,z values
        # from the coordinate data.

        MINMAX_BOUNDS=`awk 'BEGIN{first=1}{if (NF>=3){x=$1; y=$2; z=$3; if (first) {first=0; xmin=x; xmax=x; ymin=y; ymax=y; zmin=z; zmax=z;} else {if (x<xmin) xmin=x; if (x>xmax) xmax=x; if (y<ymin) ymin=y; if (y>ymax) ymax=y; if (z<zmin) zmin=z; if (z>zmax) zmax=z;}}} END{print xmin" "xmax" "ymin" "ymax" "zmin" "zmax;}' < "$tmp_atom_coords"`

        # ...and add a narrow margin (10%) around the boundaries:
        BOXSIZE_MINX=`echo $MINMAX_BOUNDS | awk '{margin=0.1; width=$2-$1; print $1-0.5*margin*width}'`
        BOXSIZE_MAXX=`echo $MINMAX_BOUNDS | awk '{margin=0.1; width=$2-$1; print $2+0.5*margin*width}'`
        BOXSIZE_MINY=`echo $MINMAX_BOUNDS | awk '{margin=0.1; width=$4-$3; print $3-0.5*margin*width}'`
        BOXSIZE_MAXY=`echo $MINMAX_BOUNDS | awk '{margin=0.1; width=$4-$3; print $4+0.5*margin*width}'`
        BOXSIZE_MINZ=`echo $MINMAX_BOUNDS | awk '{margin=0.1; width=$6-$5; print $5-0.5*margin*width}'`
        BOXSIZE_MAXZ=`echo $MINMAX_BOUNDS | awk '{margin=0.1; width=$6-$5; print $6+0.5*margin*width}'`
    else
        # By default, choose some reasonably large box:
        BOXSIZE_MINX="-100.0"
        BOXSIZE_MAXX="100.0"
        BOXSIZE_MINY="-100.0"
        BOXSIZE_MAXY="100.0"
        BOXSIZE_MINZ="-100.0"
        BOXSIZE_MAXZ="100.0"
        # ...and print message scolding the user for being lazy
        echo "----------------------------------------------------------------------" >&2
        echo "---- WARNING: Unable to determine periodic boundary conditions.   ----" >&2
        echo "----           (A default cube of volume=(200.0)^3 was used.      ----" >&2
        echo "----               This is probably not what you want!)           ----" >&2
        echo "---- It is recommended that you specify your periodic boundary    ----" >&2
        echo "---- by adding a write_once(\"Data Boundary\") command to your .lt file. ----" >&2
        echo "---- For example:                                                 ----" >&2
        #echo "----------------------------------------------------------------------" >&2
        echo "----                                                              ----" >&2
        echo "----   write_once(\"Data Boundary\") {                                   ----" >&2
        echo "----     2.51  46.79 xlo xhi                                      ----" >&2
        echo "----     -4.38 35.824 ylo yhi                                     ----" >&2
        echo "----     0.3601 42.95 zlo zhi                                     ----" >&2
        echo "----   }                                                          ----" >&2
        echo "----------------------------------------------------------------------" >&2
    fi
fi




if [ -z "$TRICLINIC" ]; then
    echo "  $BOXSIZE_MINX $BOXSIZE_MAXX xlo xhi" >> "$OUT_FILE_DATA"
    echo "  $BOXSIZE_MINY $BOXSIZE_MAXY ylo yhi" >> "$OUT_FILE_DATA"
    echo "  $BOXSIZE_MINZ $BOXSIZE_MAXZ zlo zhi" >> "$OUT_FILE_DATA"
else
    echo "triclinic parameters: XY XZ YZ = $BOXSIZE_XY $BOXSIZE_XZ $BOXSIZE_YZ" >&2
    echo "" >&2
    # Otherwise, this is a triclinic (non orthoganal) crystal basis.
    # LAMMPS represents triclinic symmetry using a different set of parameters
    # (lx,ly,lz,xy,xz,yz) than the PDB file format (alpha,beta,gamma).
    echo "  $BOXSIZE_MINX $BOXSIZE_MAXX xlo xhi" >> "$OUT_FILE_DATA"
    echo "  $BOXSIZE_MINY $BOXSIZE_MAXY ylo yhi" >> "$OUT_FILE_DATA"
    echo "  $BOXSIZE_MINZ $BOXSIZE_MAXZ zlo zhi" >> "$OUT_FILE_DATA"
    #echo "  0.000000 $BOXSIZE_X xlo xhi" >> "$OUT_FILE_DATA"
    #echo "  0.000000 $BOXSIZE_Y ylo yhi" >> "$OUT_FILE_DATA"
    #echo "  0.000000 $BOXSIZE_Z zlo zhi" >> "$OUT_FILE_DATA"
    echo "  $BOXSIZE_XY $BOXSIZE_XZ $BOXSIZE_YZ xy xz yz" >> "$OUT_FILE_DATA"
    #echo "a,b,c,alpha,beta,gamma = $BOXSIZE_A,$BOXSIZE_B,$BOXSIZE_C,$ALPHA,$BETA,$GAMMA"
fi
echo "" >> "$OUT_FILE_DATA"





if [ -s "$data_masses" ]; then
    echo "Masses" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_masses" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
else
    echo "WARNING: missing file \"$data_masses\"" >&2
fi

if [ -s "$data_pair_coeffs" ]; then
    echo "Pair Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_pair_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    PAIR_COEFFS_IN_DATA="true"
fi
if [ -s "$data_pairij_coeffs" ]; then
    echo "PairIJ Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_pairij_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    PAIR_COEFFS_IN_DATA="true"
fi
if [ -n "$PAIR_COEFFS_IN_DATA" ]; then
    if [ ! -s "$in_settings" ] || (! grep -q pair_coeff "$in_settings"); then
        echo "WARNING: no pair coeffs have been set!" >&2
    fi
fi



if [ -s "$data_bond_coeffs" ]; then
    echo "Bond Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_bond_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
else
    if [ -n "$NBONDTYPES" ] && ( [ ! -s "$in_settings" ] || (! grep -q bond_coeff "$in_settings") ); then
        echo "WARNING: no bond coeff have been set!" >&2
    fi
fi


if [ -s "$data_angle_coeffs" ]; then
    echo "Angle Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_angle_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
else
    if [ -n "$NANGLETYPES" ] && ( [ ! -s "$in_settings" ] || (! grep -q angle_coeff "$in_settings") ); then
        echo "WARNING: no angle coeffs have been set!" >&2
    fi
fi

if [ -s "$data_dihedral_coeffs" ]; then
    echo "Dihedral Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_dihedral_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
else
    if [ -n "$NDIHEDRALTYPES" ] && ( [ ! -s "$in_settings" ] || (! grep -q dihedral_coeff "$in_settings") ); then
        echo "WARNING: no dihedral coeffs have been set!" >&2
    fi
fi

if [ -s "$data_improper_coeffs" ]; then
    echo "Improper Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_improper_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
else
    if [ -n "$NIMPROPERTYPES" ] && ( [ ! -s "$in_settings" ] || (! grep -q improper_coeff "$in_settings") ); then
        echo "WARNING: no improper coeffs have been set!" >&2
    fi
fi


# data file sections specific to class2 force-fields:

if [ -s "$data_bondbond_coeffs" ]; then
    echo "BondBond Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_bondbond_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
fi

if [ -s "$data_bondangle_coeffs" ]; then
    echo "BondAngle Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_bondangle_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
fi

if [ -s "$data_middlebondtorsion_coeffs" ]; then
    echo "MiddleBondTorsion Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_middlebondtorsion_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
fi

if [ -s "$data_endbondtorsion_coeffs" ]; then
    echo "EndBondTorsion Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_endbondtorsion_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
fi

if [ -s "$data_angletorsion_coeffs" ]; then
    echo "AngleTorsion Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_angletorsion_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
fi

if [ -s "$data_angleangletorsion_coeffs" ]; then
    echo "AngleAngleTorsion Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_angleangletorsion_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
fi

if [ -s "$data_bondbond13_coeffs" ]; then
    echo "BondBond13 Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_bondbond13_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
fi

if [ -s "$data_angleangle_coeffs" ]; then
    echo "AngleAngle Coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_angleangle_coeffs" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
fi


if [ -s "$data_atoms" ]; then
    if [ -n "$ATOM_STYLE" ]; then
        echo "Atoms  # $ATOM_STYLE" >> "$OUT_FILE_DATA"
    else
        echo "Atoms  # full" >> "$OUT_FILE_DATA"
    fi
    #if [ -s "$tmp_atom_coords" ]; then
    #    echo "# (Note: x,y,z coordinates may overlap and can be modified later.)" >> "$OUT_FILE_DATA"
    #else
    #    echo "" >> "$OUT_FILE_DATA"
    #fi
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_atoms" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
else
    echo "WARNING: missing file \"$data_atoms\"" >&2
fi

if [ -s "$data_ellipsoids" ]; then
    echo "Ellipsoids" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"

    # Replace the Quaternion section, if a dump has been used to write the DATA file.
    if [ -s "$tmp_ellips_quat" ]; then 

       NATOMS=`awk 'END{print NR}' "$data_ellipsoids"`
       NATOMCRDS=`awk 'END{print NR}' "$tmp_ellips_quat"`
       if [ $NATOMS -ne $NATOMCRDS ]; then
           echo "Error: Number of atoms in coordinate file provided by user ($NATOMCRDS)" >&2
           echo "does not match the number of atoms generated in ttree file ($NATOMS)" >&2
           exit 14
       fi

       # Merge the two files.
       paste -d " " "$data_ellipsoids" "$tmp_ellips_quat" |\
       awk '{print $1" "$2" "$3" "$4" "$9" "$10" "$11" "$12}' >> "$OUT_FILE_DATA"
       echo "" >> "$OUT_FILE_DATA"

    else
       cat "$data_ellipsoids" >> "$OUT_FILE_DATA"
       echo "" >> "$OUT_FILE_DATA"
    fi
fi

if [ -s "$data_triangles" ]; then
    echo "Triangles" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_triangles" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
fi

if [ -s "$data_lines" ]; then
    echo "Lines" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_lines" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
fi

if [ -s "$data_velocities" ]; then
    echo "Velocities" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_velocities" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
#else
#    echo "Velocities" >> "$OUT_FILE_DATA"
#    echo "" >> "$OUT_FILE_DATA"
#    awk '{print $1 " 0.0 0.0 0.0"}' < "$data_atoms" >> "$OUT_FILE_DATA"
#    echo "" >> "$OUT_FILE_DATA"
fi

if [ -s "$data_bonds" ]; then
    echo "Bonds" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_bonds" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
#else
#    echo "WARNING: missing file \"$data_bonds\"" >&2
fi



if [ -s "$data_angles" ]; then
    echo "Angles" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_angles" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
#else
#    echo "WARNING: missing file \"$data_angles\"" >&2
fi

if [ -s "$data_dihedrals" ]; then
    echo "Dihedrals" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_dihedrals" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
#else
#    echo "WARNING: missing file \"$data_dihedrals\"" >&2
fi

if [ -s "$data_impropers" ]; then
    echo "Impropers" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_impropers" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
#else
#    echo "WARNING: missing file \"$data_impropers\"" >&2
fi

# OMR: if Moltemplate has been invoked without creating any atom,
# delete the spurious data file before proceeding.
if [[ $NATOMTYPES -eq 0 ]]; then
  rm -f $OUT_FILE_DATA
fi


rm -f $OUT_FILE_EXAMPLE_SCRIPT

if [ -s "$in_init" ]; then
    echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
    #  Remove leading blank spaces in the post-processed script.
    sed 's/^ *//; s/^	*//' "$in_init" > $OUT_FILE_INIT
    echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# ----------------- Init Section -----------------" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "include \"$OUT_FILE_INIT\"" >> $OUT_FILE_EXAMPLE_SCRIPT
    #echo "# \"$in_init\" typically contains various styles, dimensions, and units:" >> $OUT_FILE_EXAMPLE_SCRIPT
    #echo "include \"$in_init\"" >> $OUT_FILE_EXAMPLE_SCRIPT
    #cat "$in_init" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
fi

if [[ -e $OUT_FILE_DATA ]]; then
    echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# ----------------- Atom Definition Section -----------------" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "read_data \"$OUT_FILE_DATA\"" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
fi


if [ -s "$in_settings" ]; then
  #echo "# \"$in_settings\" typically contains coeffs, fixes, groups & modify commands:" >> $OUT_FILE_EXAMPLE_SCRIPT
  #echo "include \"$in_settings\"" >> $OUT_FILE_EXAMPLE_SCRIPT
  #cat "$in_settings" >> $OUT_FILE_EXAMPLE_SCRIPT
  if [ -z $SETTINGS_MOLC ]; then
    # Commenting out the next line.  We don't want to copy over the old file
    #cp -f "$in_settings" $OUT_FILE_SETTINGS
    # Instead we print out a warning message:
    if [ -s ${OUT_FILE_SETTINGS} ]; then
      echo "WARNING: Your files contain \"write()\" or \"write_once()\" commands which" >&2
      echo "         explicitly write text to two different but equivalent files:" >&2
      echo "           \"$in_settings\" and \"${OUT_FILE_SETTINGS}\"" >&2
      echo "         After post-processing, all of the text in either file will be copied" >&2
      echo "         into the \"${OUT_FILE_SETTINGS}\" file." >&2
      echo "         However, the order of lines in the file may be incorrect." >&2
      echo "         You should choose one file name or the other, but not both." >&2
    fi
    #...and append the contents to the end of the target file.
    #  Also, remove leading blank spaces.
    sed 's/^ *//; s/^	*//' "$in_settings" >> "${OUT_FILE_SETTINGS}"
  else
    molc.sh "$in_settings" "$in_init" >> $OUT_FILE_SETTINGS
    export MOLTEMPLATE_CITE_LIST=`printf "$MOLTEMPLATE_CITE_LIST\nRicci et al. Phys.Chem.Chem.Phys 2019 (https://doi.org/10.1039/c9cp04120f)\n"`
  fi
  echo "# ----------------- Settings Section -----------------" >> $OUT_FILE_EXAMPLE_SCRIPT
  echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
  echo "include \"$OUT_FILE_SETTINGS\"" >> $OUT_FILE_EXAMPLE_SCRIPT
  echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
fi


if [ -s "$tmp_atom_coords" ]; then
    NATOMS=`awk 'BEGIN{n=0} /^\\\$\/atom:/{n++}END{print n}' < ttree_assignments.txt`
    NATOMS_SP=`awk 'BEGIN{n=0} /^\\\${\/atom:/{n++}END{print n}' < ttree_assignments.txt`
    NATOMS=$((NATOMS + NATOMS_SP))
    NATOMCRDS=`awk '{if (NF>=3) natom+=1} END{print(natom)}' < "$tmp_atom_coords"`
    if [ $NATOMS -ne $NATOMCRDS ]; then
        echo "Error: Number of atoms in coordinate file provided by user ($NATOMCRDS)" >&2
        echo "does not match the number of atoms generated in ttree file ($NATOMS)" >&2
        exit 14
    fi

    # Copy the coordinates in $tmp_atom_coords into $OUT_FILE_DATA
    rm -f "$OUT_FILE_COORDS"
    if ! eval $PYTHON_COMMAND "${PY_SCR_DIR}/raw2data.py -ignore-atom-id " $ATOM_STYLE_ARG "$OUT_FILE_DATA" < "$tmp_atom_coords" > "$OUT_FILE_COORDS"; then
        ERR_INTERNAL
    fi
    mv -f "$OUT_FILE_COORDS" "$OUT_FILE_DATA"
    echo "copied atomic coordinates into $OUT_FILE_DATA"

    # Previously, in earlier versions of moltemplate, we used to
    # create a new input script containing "set" commands which the
    # user was supposed to tell LAMMPS to read using an "include" command.
    # If for some reason, you want to go back to doing it that way, then
    # uncomment the following 6 lines:
    #
    #rm -f "$OUT_FILE_COORDS"
    #awk '{if (NF>=3) {natom++; print "set atom "natom"  x "$1"  y "$2"  z "$3"  image 0 0 0"}}' < "$tmp_atom_coords" >> "$OUT_FILE_COORDS"
    #echo "# Load the atom coordinates:" >> $OUT_FILE_EXAMPLE_SCRIPT
    #echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
    #echo "include \"$OUT_FILE_COORDS\"" >> $OUT_FILE_EXAMPLE_SCRIPT
    #echo "" >> $OUT_FILE_EXAMPLE_SCRIPT

else
    rm -f "$OUT_FILE_COORDS"
    #echo "Warning: (moltemplate.sh)" >&2
    #echo "         Atomic coordinates were not supplied externally" >&2
    #echo "         (for example using the \"-pdb\" or \"-xyz\" arguments)." >&2
    #echo "         Hopefully you used rot(), trans() lttree commands to move" >&2
    #echo "         molecules to non-overlapping positions." >&2
fi


# Later on, it will be useful to keep track of the files that the user created.
# This must be inferred from the files ending in ".template".
# (Do this now before we discard all of the ".template" files.)

RENDERED_FILES=""
for f in *.template; do

    # Moltemplate (ttree.py) generates two copies of each file that the
    # user creates using the "write" and "write_once" commands:
    # 1) The original file requested by the user.  This file is a "rendered"
    #    file.  All of the $-style and @-style variables have been replaced
    #    with other text (usually integers).
    # 2) A version of that file that contains the original variable names
    #    instead of integers.  (These file names end in ".template".)
    # We no longer care about the files ending in ".template", but we can loop
    # over files ending in ".template" created by moltemplate.  We can infer
    # the name of the other file name by removing the ".template" suffix.

    bn=`basename "$f" .template`  #file name without the ".template" suffix
    if [ ! -f "$bn" ]; then
        continue
    fi
    if [ -n "$RENDERED_FILES" ]; then
        printf -v RENDERED_FILES "${RENDERED_FILES}\n$bn"
    else
        RENDERED_FILES="$bn"
    fi
done



# ############## CLEAN UP ################

# A lot of files have been created along the way.
# However only a few of them are actually useful.
#
# Optional: clean up some non-essential temporary
#           files created by running ttree
# We move the non-essential files into a different directory
# (but we don't delete them).

if [ ! -d output_ttree ]; then
    mkdir output_ttree
fi




# Move temporary files into the "output_ttree/" directory:
IFS=$CR
rm -f ttree_replacements.txt >/dev/null 2>&1 || true
for file in $MOLTEMPLATE_TEMP_FILES; do
    if [ -e "$file" ]; then
        rm -f "output_ttree/$file" >/dev/null 2>&1 || true
        #echo "file=\"$file\""
        mv "$file" output_ttree/ >/dev/null 2>&1 || true
    fi
done
IFS=$OIFS


# We no longer need the "ttree_assignments_static.txt" file.
rm -f ttree_assignments_static.txt



# ############## DEAL WITH CUSTOM NON-STANDARD SECTIONS ################


# N_data_prefix=`expr length "$data_prefix"` <-- not posix compliant. AVOID
N_data_prefix=${#data_prefix}  #<-- works even if $data_prefix contains spaces
#for file_name in "${data_prefix}"*; do
ls "${data_prefix}"* 2> /dev/null | while read file_name; do
    #If using bash:
    #SECTION_NAME="${file_name:$N_data_prefix}"
    #If using sh:
    #SECTION_NAME=`expr substr "$file_name" $(($N_data_prefix+1)) 1000000` <-- not posix compliant. AVOID
    SECTION_NAME=`echo "" | awk -v FN="$file_name" -v N="$N_data_prefix" 'END{print substr(FN,N+1,1000000)}'`
    # Create a new section in the data file
    # matching the portion of the name of
    # the file after the data_prefix.
    echo "" >> "$OUT_FILE_DATA"
    echo "$SECTION_NAME" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    cat "$file_name" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    mv -f "$file_name" output_ttree/
done


if [ -e "$data_prefix_no_space" ]; then
    echo "" >> "$OUT_FILE_DATA"
    cat "$data_prefix_no_space" >> "$OUT_FILE_DATA"
    echo "" >> "$OUT_FILE_DATA"
    mv -f "$data_prefix_no_space" output_ttree/
fi



if [ -e "$OUT_FILE_DATA" ]; then
    mv -f "$OUT_FILE_DATA" "$OUT_FILE_DATA.tmp" >/dev/null 2>&1 || true
    #dos2unix < "$OUT_FILE_DATA.tmp" > "$OUT_FILE_DATA"
    tr -d '\r' < "$OUT_FILE_DATA.tmp" > "$OUT_FILE_DATA"
    rm -f "$OUT_FILE_DATA.tmp" >/dev/null 2>&1 || true
fi



# Convert custom files with names like "In Charges" to "system.in.charges"


#N_in_prefix=`expr length "$in_prefix"` <-- not posix compliant. AVOID.
N_in_prefix=${#in_prefix}  #<-- works even if $in_prefix contains spaces
#for file_name in "${in_prefix}"*; do
ls "${in_prefix}"* 2> /dev/null | while read file_name; do
    #If using bash:
    #SECTION_NAME="${file_name:$N_in_prefix}"
    #If using sh:
    #SECTION_NAME=`expr substr "$file_name" $(($N_in_prefix+1)) 1000000` <-- not posix compliant. AVOID
    SECTION_NAME=`echo "" | awk -v FN="$file_name" -v N="$N_in_prefix" 'END{print substr(FN,N+1,1000000)}'`
    #FILE_SUFFIX=`echo "$SECTION_NAME" | awk '{print tolower($0)}'`
    FILE_SUFFIX=`echo "$SECTION_NAME" | awk '{print tolower($0)}'|sed 's/ /./g'`
    # Create a new section in the lammps input script
    # matching the portion of the name of
    # the file after the in_prefix.
    echo "" >> "$OUT_FILE_EXAMPLE_SCRIPT"
    echo "# ----------------- $SECTION_NAME Section -----------------" >> $OUT_FILE_EXAMPLE_SCRIPT

    # Commenting out the next line.
    #
    #cp -f "$file_name" "${OUT_FILE_EXAMPLE_SCRIPT}.${FILE_SUFFIX}"
    #
    # We don't want overwrite "${OUT_FILE_EXAMPLE_SCRIPT}.${FILE_SUFFIX}"
    # with the contents of "$file_name" because it's possible that the user
    # might have created both files (eg "In Charges" and "system.in.charges"),
    # and we don't want to erase anything.  Instead, print a warning message
    if [ -s "${OUT_FILE_INPUT_SCRIPT}.${FILE_SUFFIX}" ]; then
      echo "WARNING: Your files contain \"write()\" or \"write_once()\" commands which" >&2
      echo "         explicitly write text to two different but equivalent files:" >&2
      echo "           \"$file_name\" and \"${OUT_FILE_INPUT_SCRIPT}.${FILE_SUFFIX}\"" >&2
      echo "         After post-processing, all of the text in either file will be copied" >&2
      echo "         into the \"${OUT_FILE_INPUT_SCRIPT}.${FILE_SUFFIX}\" file." >&2
      echo "         However, the order of lines in the file may be incorrect." >&2
      echo "         You should choose one file name or the other, but not both." >&2
    fi
    #...and append the contents to the end of the target file.
    #  Also, remove leading blank spaces.
    sed 's/^ *//; s/^	*//' "$file_name" >> "${OUT_FILE_INPUT_SCRIPT}.${FILE_SUFFIX}"
    echo "" >> "$OUT_FILE_EXAMPLE_SCRIPT"
    echo "include \"${OUT_FILE_INPUT_SCRIPT}.${FILE_SUFFIX}\"" >> $OUT_FILE_EXAMPLE_SCRIPT
    mv -f "$file_name" output_ttree/
done

if [ -e "$in_prefix_no_space" ]; then
    echo "" >> "$OUT_FILE_EXAMPLE_SCRIPT"
    cat "$in_prefix_no_space" >> "$OUT_FILE_EXAMPLE_SCRIPT"
    echo "" >> "$OUT_FILE_EXAMPLE_SCRIPT"
    mv -f "$in_prefix_no_space" output_ttree/
fi

# ------ Postprocess coeff commands after rendering. ----------
# Swap the order of atom types I, J in all "pair_coeff I J ..." commands
# whenever I > J.  Do this for every input script file generated by moltemplate.
# (Perhaps later I'll check to make sure the user did not specify contradictory
#  "pair_coeff I J ..." and "pair_coeff J I ..." commands, but I won't do it
#  here, because at this point we've thrown away the original atom type names,
#  so there's no easy way to explain the problem to the user if there is one.)

# -- First we must figure out which files we need to consider. --
# Ugly details.  We already created a variable ($RENDERED_FILES) storing a
# list of files which are created by the user.  However some of those files no
# longer exist. We also need do add several files to this list and remove
# duplicates.  Lacking data structures, SH and BASH are really bad at this kind
# thing, so the code is ugly.  Must create temporary files and sort them.

echo "$RENDERED_FILES" > RENDERED_FILENAMES.tmp
echo "$OUT_FILE_INIT" >> RENDERED_FILENAMES.tmp
echo "$OUT_FILE_EXAMPLE_SCRIPT " >> RENDERED_FILENAMES.tmp
echo "$OUT_FILE_SETTINGS " >> RENDERED_FILENAMES.tmp
echo "$RENDERED_FILES_WITH_COEFF_COMMANDS" >> RENDERED_FILENAMES.tmp

# Now remove duplicate file names
sort RENDERED_FILENAMES.tmp > RENDERED_FILENAMES.tmp2
awk '{if ((NF>0)&&(NR>1)&&($0!=prev)) {print $0} prev=$0}' \
    < RENDERED_FILENAMES.tmp2 > RENDERED_FILENAMES.tmp3

# Now update the list of files stored in $RENDERED_FILES
RENDERED_FILES=""
while read file_name; do
    if [ ! -f "$file_name" ]; then
        continue  # make sure the file exists
    fi
    if [ -n "$RENDERED_FILES" ]; then
        printf -v RENDERED_FILES "${RENDERED_FILES}\n$file_name"
    else
        RENDERED_FILES="$file_name"
    fi
done < RENDERED_FILENAMES.tmp3

rm -f RENDERED_FILENAMES.tmp*


# Now filter out files which do not contain "_coeff" commands
RENDERED_FILES_WITH_COEFF_COMMANDS=""
IFS=$CR
for file_name in $RENDERED_FILES; do
    if [ ! -f "$file_name" ]; then
        continue  # make sure the file exists
    fi
    HAS_COEFF_COMMANDS=`awk 'BEGIN{has_coeffs=0} {if ((NF>0)&&(($1=="pair_coeff")||($1=="bond_coeff")||($1=="angle_coeff")||($1=="dihedral_coeff")||($1=="improper_coeff"))) has_coeffs=1} END{print has_coeffs}' < "$file_name"`
    if [ "$HAS_COEFF_COMMANDS" -eq "1" ]; then
        if [ -n "$RENDERED_FILES_WITH_COEFF_COMMANDS" ]; then
            printf -v RENDERED_FILES_WITH_COEFF_COMMANDS "${RENDERED_FILES_WITH_COEFF_COMMANDS}\n$file_name"
        else
            RENDERED_FILES_WITH_COEFF_COMMANDS="$file_name"
        fi
    fi
done
IFS=$OIFS



# ---- Now loop over all user-created files containing "coeff" commands. ----
echo "" > input_scripts_so_far.tmp
IFS=$CR
for file_name in $RENDERED_FILES_WITH_COEFF_COMMANDS; do
    if [ -s "$file_name" ]; then
        echo "postprocessing file \"$file_name\"" >&2
        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/postprocess_input_script.py" input_scripts_so_far.tmp < "$file_name" > "$file_name.tmp"; then
            ERR_INTERNAL
        fi
        echo "" >&2

        mv -f "$file_name.tmp" "$file_name"
        #cat "$file_name" >> input_scripts_so_far.tmp
        #dos2unix < "$file_name" >> input_scripts_so_far.tmp
        tr -d '\r' < "$file_name" >> input_scripts_so_far.tmp

        # Delete all "bond_style" statements when no bond types are defined
        if [ -z "$NBONDTYPES" ]; then
            awk '{if ($1!="bond_style") print $0}' < "$file_name" > "${file_name}.tmp"
            mv -f "$file_name.tmp" "$file_name"
        fi
        # Delete all "angle_style" statements when no angle types are defined
        if [ -z "$NANGLETYPES" ]; then
            awk '{if ($1!="angle_style") print $0}' < "$file_name" > "${file_name}.tmp"
            mv -f "$file_name.tmp" "$file_name"
        fi
        # Delete all "dihedral_style" statements when no dihedral types are defined
        if [ -z "$NDIHEDRALTYPES" ]; then
            awk '{if ($1!="dihedral_style") print $0}' < "$file_name" > "${file_name}.tmp"
            mv -f "$file_name.tmp" "$file_name"
        fi
        # Delete all "improper_style" statements when no improper types are defined
        if [ -z "$NIMPROPERTYPES" ]; then
            awk '{if ($1!="improper_style") print $0}' < "$file_name" > "${file_name}.tmp"
            mv -f "$file_name.tmp" "$file_name"
        fi
    fi
done
IFS=$OIFS


# Process any custom input script files created by the user afterwards
#ls "${OUT_FILE_INPUT_SCRIPT}"* 2> /dev/null | while read file_name; do
for file_name in "${OUT_FILE_INPUT_SCRIPT}."*; do
    if [ ! -f "$file_name" ]; then
        # Special case: If there aren't any files which match the pattern
        # ("${in_prefix}"*.template), then the bash glob expansion will 
        # generate a list with only one entry ("${in_prefix}"*.template).
        # In that case, there will be no file of that name, and we should skip
        # this entry (or break out of the loop, since there's only one entry)
        continue
    fi
    if [ "$file_name" = "$OUT_FILE_INIT" ] || [ $file_name = "$OUT_FILE_INPUT_SCRIPT" ] || [ $file_name = "$OUT_FILE_SETTINGS" ]; then
        continue
    fi
    if [[ "$file_name" == *.tmp ]]; then
        continue
    fi
    if [[ "$file_name" == *.template ]]; then
        continue
    fi
    echo "postprocessing file \"$file_name\"" >&2
    if ! $PYTHON_COMMAND "${PY_SCR_DIR}/postprocess_input_script.py" input_scripts_so_far.tmp < "$file_name" > "$file_name.tmp"; then
        ERR_INTERNAL
    fi
    echo "" >&2
    mv "$file_name.tmp" "$file_name"
    cat "$file_name" >> input_scripts_so_far.tmp
done

rm -f input_scripts_so_far.tmp




# ############ Optional: Add a fake run section as an example ############
if [ ! -z $APPEND_EXAMPLE_SCRIPT ]; then
    echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# ----------------- Run Section -----------------" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# The lines above define the system you want to simulate." >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# What you do next is up to you." >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# Typically a user would minimize and equilibrate" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# the system using commands similar to the following:" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "#  ----   examples   ----" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "#" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "#  -- minimize --" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# minimize 1.0e-5 1.0e-7 1000 10000" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# (Note: Some fixes, for example \"shake\", interfere with the minimize command," >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "#        You can use the \"unfix\" command to disable them before minimization.)" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "#  -- declare time step for normal MD --" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# timestep 1.0" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "#  -- run at constant pressure (Nose-Hoover)--" >> $OUT_FILE_EXAMPLE_SCRIPT
    #echo "# timestep 1.0" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# fix   fxnpt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 1.0">>$OUT_FILE_EXAMPLE_SCRIPT
    echo "#  -- ALTERNATELY, run at constant volume (Nose-Hoover) --" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# fix   fxnvt all nvt temp 300.0 300.0 500.0 tchain 1" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "#  -- ALTERNATELY, run at constant volume using Langevin dynamics. --" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "#  -- (This is good for sparse CG polymers in implicit solvent.)   --" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# fix fxLAN all langevin 300.0 300.0 5000 48279" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# fix fxNVE all nve  #(<--needed by fix langevin)" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "#  -- Now, finally run the simulation --" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "# run   50000" >> $OUT_FILE_EXAMPLE_SCRIPT
    #echo "# write_restart system_after_nvt.rst" >> $OUT_FILE_EXAMPLE_SCRIPT
    #echo "# run   50000" >> $OUT_FILE_EXAMPLE_SCRIPT
    #echo "# write_restart system_after_npt.rst" >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "#  ---- (end of examples) ----">> $OUT_FILE_EXAMPLE_SCRIPT
    #echo "# It is the responsibility of the user to learn LAMMPS and specify these">>$OUT_FILE_EXAMPLE_SCRIPT
    #echo "# these commands." >> $OUT_FILE_EXAMPLE_SCRIPT
    echo "" >> $OUT_FILE_EXAMPLE_SCRIPT
fi


# Finally, if the -vmd argument was included, start up VMD and
# view the system (using topotools to convert the

if [ ! -z $RUN_VMD_AT_END ]; then

    echo "topo readlammpsdata $OUT_FILE_DATA $ATOM_STYLE" > vmd_viz_moltemplate.tcl.tmp
    bn=`basename $OUT_FILE_DATA .data`
    echo "animate write psf $bn.psf" >> vmd_viz_moltemplate.tcl.tmp
    vmd -e vmd_viz_moltemplate.tcl.tmp
    rm -rf vmd_viz_moltemplate.tcl.tmp

fi


# Lastly, ask users to cite our paper.
# Hopefully this message is not too intrusive.
echo   "-------------------------------------------------------------" >&2
printf "If this software is useful in your research, please cite\n" >&2
if [ `printf "$MOLTEMPLATE_CITE_LIST\n" | awk 'END{print NR}'` -gt 1 ]; then
  i_CITE=1
  printf "$MOLTEMPLATE_CITE_LIST\n" | while read CITATION; do
    echo "${i_CITE}) $CITATION" >&2
    i_CITE=$(( i_CITE+1 ))
  done
else
  printf "$MOLTEMPLATE_CITE_LIST\n" >&2
fi
echo   "-------------------------------------------------------------" >&2
