crysgen

#!/usr/bin/env bash
# Brief: Generates a beta-pvdf crystal from a given periodic n-mer (maybe energy minimized output of polygen).
## User input required: A properly aligned single periodic n-mer, na (number of molecules along a), nb, nc and b(nm)(unit cell dimension along dipole). 
# Ref: Zhu et al. 2008, Comput. Mat. Sci. 44, 224-229
# Usage: crysgen <path/to/.gro> <nx> <ny> <nz> [<b>]  # b default 0.491nm

set -e
tmpdir="$(mktemp -u)"; mkdir "${tmpdir}"

# Warn user
echo "Reminder: Is your molecule properly aligned? Dipoles along Y and long axis along Z? (Hint: Use output of polygen)" >&2

# Inputs
mol_gro="${1}" # .gro file of single n-mer beta pvdf
[[ -n "${mol_gro}" ]] && [[ -f "${mol_gro}" ]] || { echo "Provide .gro of single properly aligned n-mer"; exit 1;}
res_per_mol="$(($(sed -n 2p "${mol_gro}")/6))" # Number of residues per molecule from <#atoms_per_molecule>/<#atoms_per_residue>
na="${2}"; nb="${3}"; nc="${4}" # Number of molecules along a, b and c directions. Total number 2 x na x nb x nc (Twice because one + one staggered along a,b)
( [[ "${na}" ]] && [[ "${nb}" ]] && [[ "${nc}" ]] ) || { echo "Provide number of mols along a,b,c"; exit 1;}
b="${5:-0.491}" # b in nm. Give all the precision you've got. Default : 0.491nm.

# Parameters:
declare -r a=0.850 c=0.256
declare -r a_half="$(echo "${a}/2" | bc -l)" b_half="$(echo "${b}/2" | bc -l)" # Might also use "$(echo "scale=5;${b}/2" | bc)"
declare -r box_a="$(echo "${na}*${a}" | bc)" box_b="$(echo "${nb}*${b}" | bc)"
declare -r box_c="$(echo "${nc}*${res_per_mol}*${c}" | bc)"

# Build original grid from single molecules
sed '$c0.0 0.0 0.0' "${mol_gro}" > "${tmpdir}"/nobox_single_mol.gro
gmx genconf -norot -nbox "${na}" "${nb}" "${nc}" -dist "${a}" "${b}" "${c}" -f "${tmpdir}"/nobox_single_mol.gro -o "${tmpdir}"/orig_grid.gro

# Shift the original grid in the X-Y or a-b plane to get the staggered grid
gmx editconf -translate "${a_half}" "${b_half}" 0 -f "${tmpdir}"/orig_grid.gro -o "${tmpdir}"/staggered_grid.gro 

# Superimpose both grids to get the quasi hexagonal lattice
natoms_orig="$(sed -n 2p "${tmpdir}"/orig_grid.gro)"
natoms="$((natoms_orig*2))"
echo -e "PVDF\n${natoms}" > "${tmpdir}"/super_grid.gro
head -n -1 "${tmpdir}"/orig_grid.gro | tail -n +3 >> "${tmpdir}"/super_grid.gro
head -n -1 "${tmpdir}"/staggered_grid.gro | tail -n +3 >> "${tmpdir}"/super_grid.gro
echo "${box_a} ${box_b} ${box_c}" >> "${tmpdir}"/super_grid.gro # Perfect box size to make the crystal
# Centre the whole system within box
gmx editconf -f "${tmpdir}"/super_grid.gro -o crystal.gro -bt triclinic -angles 90.0 90.0 90.0 -box "${box_a}" "${box_b}" "${box_c}" -resnr 1
# Note the -resnr option above, to renumber the residues

# Cleanup
rm -rf "${tmpdir}"