#!/bin/bash
#########################################################
# Script to perform a geometry optimistion using lammps #
#########################################################

set -e

#########################
# Test usage is correct #
#########################

if [ $# -ne 3 ]; then
         echo 1>&2 "Usage: lammps_relax <exe> <pressure> <seed>"
         exit 127
fi

###########################
# Set the input variables #
###########################

exe=$1
setextpress=$2
seed=$3

root=${seed%%-*}
root=${root##*\/}

# Build the lammps input file

cat << _EOF > $seed.in
#  --------------------- Initialisation  ---------------------
units real # distances in A, energy in Kcal/mole, pressure in atmospheres
dimension 3
boundary p p p          # periodic in x, y, z
atom_style   charge

#  --------------------- set pressure  ---------------------
variable extpress equal 9869.2327*$setextpress # Convert GPa to atm

#  --------------------- Atom definition  ---------------------
box tilt large          # relax restictions on box tilt factors
read_data $seed.conf

neigh_modify every 1 delay 0 check yes # update neigbour every step
neighbor 2.0 bin                       # neighbour list skin

#  --------------------- Potential definition  ---------------------
include $seed.pp

#  --------------------- record traj  ---------------------
dump mintraj all atom 1 $seed.lammpstrj          # dump the trajectory
dump_modify mintraj sort id                 # sort atoms by numerical label

#   --------------------- initial Relaxation of the structure   --------------------- 
min_style cg
min_modify line quadratic

# Multiple minimisations to prevent trapping
fix minbox all box/relax tri \${extpress} vmax 0.1
minimize 1.0e-8 1.0e-6 10000 1000000    
fix minbox all box/relax tri \${extpress} vmax 0.01
minimize 1.0e-10 1.0e-6 10000 1000000
fix minbox all box/relax tri \${extpress} vmax 0.001
minimize 1.0e-12 1.0e-8 10000 100000
fix minbox all box/relax tri \${extpress} vmax 0.0001
minimize 1.0e-12 1.0e-8 10000 100000
unfix minbox

#   ---------------------  collect results   --------------------- 
# Variables defined for printing
variable h equal enthalpy*0.043364104            # convert kcal/mol to eV
variable e equal etotal*0.043364104              # convert kcal/mol to eV
variable v equal vol
variable p equal press*0.000101325               # convert atm back to GPa
variable a equal cella
variable b equal cellb
variable c equal cellc
variable alp equal cellalpha
variable bet equal cellbeta
variable gam equal cellgamma
# Print AIRSS-specific data in $seed.lammps
print "Lattice parameters:  \${a} \${b} \${c} \${alp} \${bet} \${gam}" file $seed.lammps
print "Volume: \$v" append $seed.lammps
print "Pressure: \$p" append $seed.lammps
print "Enthalpy: \$h" append $seed.lammps
_EOF

# Run lammps

cabal cell conf < $seed.cell > $seed.conf

lammps < $seed.in > $seed.lout

rm -f log.lammps

# Construct a fake Castep output file

grep Pressure: $seed.lout | tail -1 | awk '{print " *  Pressure: "$2}' > $seed.castep
grep Enthalpy: $seed.lout | tail -1 | awk '{print " lammps: Final Enthalpy     = "$2}' >> $seed.castep
grep Volume:   $seed.lout | tail -1 | awk '{print "Current cell volume = "$2}' >> $seed.castep

# Save the final structure in Castep -out.cell format

lammps2cell $seed > $seed-out.cell

# Return the key results

cat $seed.lout

exit 0