#!/usr/bin/env perl 

#==================================================================================#
#                                  airss.pl                                        #
#==================================================================================#
#                                                                                  #
# This file is part of the AIRSS structure prediction package.                     #
#                                                                                  #
# AIRSS is free software; you can redistribute it and/or modify it under the terms #
# of the GNU General Public License version 2 as published by the Free Software    #
# Foundation.                                                                      #
#                                                                                  #
# This program is distributed in the hope that it will be useful, but WITHOUT ANY  #
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A  #
# PARTICULAR PURPOSE.  See the GNU General Public License for more details.        #           
#                                                                                  #
# You should have received a copy of the GNU General Public License along with this#
# program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street,#                   
# Fifth Floor, Boston, MA  02110-1301, USA.                                        #
#                                                                                  #
#----------------------------------------------------------------------------------#
# This main AIRSS script, generate and run random structures generated by buildcell#
#----------------------------------------------------------------------------------#
# Written by Chris Pickard, Copyright (c) 2005-2020                                #
#----------------------------------------------------------------------------------#
#                                                                                  #
#==================================================================================#

use strict;
use Getopt::Long;
use File::Copy;
use Time::HiRes qw (time);

sub usage {
  printf STDERR "\n";
  printf STDERR "      .o.       ooooo ooooooooo.    .oooooo..o  .oooooo..o \n";
  printf STDERR "     .888.      '888' '888   'Y88. d8P'    'Y8 d8P'    'Y8 \n";
  printf STDERR "    .8:888.      888   888   .d88' Y88bo.      Y88bo.      \n";
  printf STDERR "   .8' '888.     888   888ooo88P'   ':Y8888o.   ':Y8888o.  \n";
  printf STDERR "  .88ooo8888.    888   888'88b.         ':Y88b      ':Y88b \n";
  printf STDERR " .8'     '888.   888   888  '88b.  oo     .d8P oo     .d8P \n";
  printf STDERR "o88o     o8888o o888o o888o  o888o 8::88888P'  8::88888P'  \n";
  printf STDERR "                                                           \n";
  printf STDERR "          Ab Initio Random Structure Searching             \n";
  printf STDERR "          Chris J. Pickard   (cjp20\@cam.ac.uk)            \n";
  printf STDERR "                 Copyright (c) 2005-2020                   \n";
  printf STDERR "                                                           \n";
  printf STDERR "Please cite the following:                                 \n";
  printf STDERR "                                                           \n";
  printf STDERR "[1] C.J. Pickard and R.J. Needs, PRL 97, 045504 (2006)     \n";
  printf STDERR "[2] C.J. Pickard and R.J. Needs, JPCM 23, 053201 (2011)    \n";
  printf STDERR "\n";
  printf STDERR "Contact airss\@msm.cam.ac.uk for support\n";
  printf STDERR "\n";
  printf STDERR "Usage: airss.pl [-pressure] [-prand] [-pmin] [-build] [-pp0] [-gosh] [-pp3] [-repose] [-gulp] [-lammps] [-gap] [-psi4] [-vasp] [-qe] [-python] [-castep] [-exec] [-cluster] [-slab] [-dos] [-workdir] [-max] [-num] [-amp] [-camp] [-mode] [-minmode] [-sim] [-symm] [-nosymm] [-mpinp] [-steps] [-best] [-track] [-keep] [-pack] [-harvest] [-seed]\n";
  printf STDERR "       -pressure f  Pressure (0.0)\n";
  printf STDERR "       -prand       Randomise pressures (false)\n";
  printf STDERR "       -pmin        Minimum pressure for randomisation (0.0)\n";
  printf STDERR "       -build       Build structures only (false)\n";
  printf STDERR "       -pp0         Use pair potentials rather than Castep (0D) (false)\n";
  printf STDERR "       -gosh        Use pair potentials rather than Castep (0D) (false)\n";
  printf STDERR "       -pp3         Use pair potentials rather than Castep (3D) (false)\n";
  printf STDERR "       -repose      Use data derived potentials rather than Castep (3D) (false)\n";
  printf STDERR "       -gulp        Use gulp rather than Castep (false)\n";
  printf STDERR "       -lammps      Use LAMMPS rather than Castep (false)\n";
  printf STDERR "       -gap         Use GAP through QUIP/QUIPPY/ASE (false)\n";
  printf STDERR "       -ps4         Use psi4 (false)\n";
  printf STDERR "       -vasp        Use VASP (false)\n";
  printf STDERR "       -qe          Use Quantum Espresso (false)\n";
  printf STDERR "       -python      Use relax.py script (false)\n";
  printf STDERR "       -castep      Use CASTEP in addition to alternative (false)\n";
  printf STDERR "       -exec        Use this executable\n";
  printf STDERR "       -launch      Use this parallel launcher\n";
  printf STDERR "       -cluster     Use cluster settings for symmetry finder (false)\n";
  printf STDERR "       -slab        Use slab settings (false)\n";
  printf STDERR "       -dos         Calculate DOS at Ef (false)\n";
  printf STDERR "       -workdir  s  Work directory ('.')\n";
  printf STDERR "       -max      n  Maximum number of structures (1000000)\n";
  printf STDERR "       -num      n  Number of trials (0)\n";
  printf STDERR "       -amp      f  Amplitude of ion move (-1.5)\n";
  printf STDERR "       -camp     f  Amplitude of cell move (0.0)\n";
  printf STDERR "       -mode        Choose moves based on low lying vibrational modes (false)\n";
  printf STDERR "       -minmode  n  Lowest mode (4)\n";
  printf STDERR "       -sim      f  Threshold for structure similarity (0.0)\n";
  printf STDERR "       -symm     f  Symmetrise on-the-fly (0.0)\n";
  printf STDERR "       -nosymm   f  No symmetry (0)\n";
  printf STDERR "       -mpinp    n  Number of cores per mpi Castep (1)\n";
  printf STDERR "       -steps    n  Max number of geometry optimisation steps (400)\n";
  printf STDERR "       -best        Only keep the best structures for each composition (false)\n";
  printf STDERR "       -track       Keep the track of good structures during relax and shake (false)\n";
  printf STDERR "       -keep        Keep intermediate files (false)\n";
  printf STDERR "       -pack        Concatenate the res files (false)\n";
  printf STDERR "       -harvest     Collect intermediate structures as res files (false)\n";
  printf STDERR "       -seed     s  Seedname ('NONE')\n";
  exit();
}

my ($opt_pressure,$opt_prand,$opt_pmin,$opt_build,$opt_pp0,$opt_gosh,$opt_pp3,$opt_repose,$opt_gulp,$opt_lammps,$opt_gap,$opt_psi4,$opt_vasp,$opt_qe,$opt_python,$opt_castep,$opt_exec,$opt_launch,$opt_cluster,$opt_slab,$opt_dos,$opt_workdir,$opt_max,$opt_num,$opt_nbest,$opt_amp,$opt_camp,$opt_mode,$opt_minmode,$opt_sim,$opt_symm,$opt_nosymm,$opt_mpinp,$opt_steps,$opt_best,$opt_track,$opt_keep,$opt_pack,$opt_harvest,$opt_seed,$opt_help)
  
  = (
    0.0,         # => $opt_pressure
    0,           # => $opt_prand
    0.0,         # => $opt_pmin
    0,           # => $opt_build
    0,           # => $opt_pp0
    0,           # => $opt_gosh
    0,           # => $opt_pp3
    0,           # => $opt_repose
    0,           # => $opt_gulp
    0,           # => $opt_lammps
    0,           # => $opt_gap
    0,           # => $opt_psi4
    0,           # => $opt_vasp
    0,           # => $opt_qe
    0,           # => $opt_python
    0,           # => $opt_castep
    "",          # => $opt_exec
    "",          # => $opt_launch
    0,           # => $opt_cluster
    0,           # => $opt_slab
    0,           # => $opt_dos
    ".",         # => $opt_workdir
    1000000,     # => $opt_max
    0,           # => $opt_num
    1,           # => $opt_nbest
    -1.5,        # => $opt_amp
    0.0,         # => $opt_camp
    0,           # => $opt_mode
    4,           # => $opt_minmode
    0.0,         # => $opt_sim
    0.0,         # => $opt_symm
    0,           # => $opt_nosymm
    1,           # => $opt_mpinp
    400,         # => $opt_steps
    0,           # => $opt_best
    0,           # => $opt_track
    0,           # => $opt_keep
    0,           # => $opt_pack
    0,           # => $opt_harvest
    "NONE",      # => $opt_seed
    0            # => $opt_help
    );

my $commandline = (join " ", @ARGV);

&GetOptions("pressure=f" => \$opt_pressure,
            "prand"      => \$opt_prand,
            "pmin=f"     => \$opt_pmin,
            "build"      => \$opt_build,
            "pp0"        => \$opt_pp0,
            "gosh"       => \$opt_gosh,
            "pp3"        => \$opt_pp3,
            "repose"     => \$opt_repose,
            "gulp"       => \$opt_gulp,
            "lammps"     => \$opt_lammps,
            "gap"        => \$opt_gap,
            "psi4"       => \$opt_psi4,
            "vasp"       => \$opt_vasp,
            "qe"         => \$opt_qe,
            "python"     => \$opt_python,
            "castep"     => \$opt_castep,
            "exec=s"     => \$opt_exec,
            "launch=s"   => \$opt_launch,
            "cluster"    => \$opt_cluster,
            "slab"       => \$opt_slab,
            "dos"        => \$opt_dos,
            "workdir=s"  => \$opt_workdir, 
            "max=n"      => \$opt_max,
            "num=n"      => \$opt_num,
            "nbest=n"    => \$opt_nbest,
            "amp=f"      => \$opt_amp,
            "camp=f"     => \$opt_camp,
            "mode"       => \$opt_mode,
            "minmode=n"  => \$opt_minmode,
            "sim=f"      => \$opt_sim,
            "symm=f"     => \$opt_symm,
            "nosymm"     => \$opt_nosymm,
            "mpinp=n"    => \$opt_mpinp,
            "steps=n"    => \$opt_steps,
            "best"       => \$opt_best,
            "track"      => \$opt_track,
            "keep"       => \$opt_keep,
            "pack"       => \$opt_pack,
            "harvest"    => \$opt_harvest,
            "seed=s"     => \$opt_seed,
            "h|help"     => \$opt_help)|| &usage();

if ($opt_help or $opt_seed eq "NONE") {
  &usage();
}

# Initialise loop counters etc

my ($high_energy,$i,$j,$k) = (999999999999999999999.999,0,0,0);

my $original_pressure = $opt_pressure;

my $best_enthalpy=$high_energy;
my $best_structure="";

my $uniquehead=$opt_seed."-".$$;
my $timestamp=time()*1000000%10000;
my $uniquenew = $uniquehead."-".$timestamp;

if ($opt_num < 0) { $opt_num=100000000000000000000000 }

# Define the executable to be used

my $executable='';

if ( $opt_vasp ) {
  $executable='vasp';
} elsif ( $opt_qe ) {
  $executable='pw.x';
} else {
  $executable='castep';
}

if ( $opt_exec ne '' ) {$executable=$opt_exec};

my $launcher='mpirun -np';

if ( $opt_launch ne '' ) {$launcher=$opt_launch};

if ( $opt_mpinp > 0 ) { $executable= '"'.$launcher.' '.$opt_mpinp.' '.$executable.'"' }

# Make work directory

system "mkdir -p $opt_workdir";

# Check that necessary files exist

if ( -e $opt_seed.".cell" ) {
  copy($opt_seed.".cell",$opt_workdir."/".$uniquenew.".cell")
} else {
  print "<seed>.cell not found\n";
  exit
}

if ( ( $opt_qe ) && !( -e $opt_seed.".qe" ) ) {
  print "QE flag used (-qe), but <seed>.qe not found\n";
  exit
}

# The main loop over random structures

if ( -e "stop" ) {unlink("stop")}

my $nt=0;

for ($i = 1; $i <= $opt_max; $i++) {

  if ( -e "stop" ) {last}
  
  $nt=$nt+1;

  # Set the pressure

  if ($opt_prand) {
    $opt_pressure=$opt_pmin+($original_pressure-$opt_pmin)*rand(1);
  } else {
    $opt_pressure=$original_pressure;
  }
  
  # Set the unique file names

  my $unique = $uniquehead."-".$timestamp."-".$i;
     
  # Make a link to the seed param or pp

  if ( $opt_pp3 || $opt_repose || $opt_lammps ) {
    copy($opt_seed.".par",$opt_workdir."/".$unique.".par");
    copy($opt_seed.".pp",$opt_workdir."/".$unique.".pp");
    copy($opt_seed.".ppp",$opt_workdir."/".$unique.".ppp");
    copy($opt_seed.".ddp",$opt_workdir."/".$unique.".ddp");
    copy($opt_seed.".eddp",$opt_workdir."/".$unique.".eddp");
  } elsif ( $opt_vasp )  {
    copy($opt_seed.".POTCAR",$opt_workdir."/".$unique.".POTCAR");
    my $vasp_pressure=10.0*$opt_pressure;
    system "(grep -v PSTRESS $opt_seed.INCAR ; echo PSTRESS = $vasp_pressure)  > $opt_workdir/$unique.INCAR";
  } elsif ( $opt_qe ) {   
    copy($opt_seed.".qe",$opt_workdir."/".$unique.".qe"); 
  } else {
    copy($opt_seed.".param",$opt_workdir."/".$unique.".param"); 
  }

  # Prepare for extra run of CASTEP
  
  if ( $opt_castep ) {
    copy($opt_seed.".param",$opt_workdir."/".$unique.".param"); 
  }
  
  # Keep a copy of the cell file

  my $md5=substr(`md5sum $opt_seed.cell`,0,8);

  copy($opt_seed.".cell",$opt_workdir."/".$opt_seed.".cell.".$md5);

  if ( ! $opt_mode || $j == 0) {
  
    # Generate a random structure based on the directives in the seed cell

    if ($opt_nbest <= 1) {

      system "buildcell < $opt_workdir/$uniquenew.cell > $opt_workdir/$unique.cell 2> /dev/null";

    } else {

      # Calculate single point energies and choose the best of n

      if ($opt_build||$opt_pp0||$opt_gosh||$opt_pp3||$opt_gulp||$opt_lammps||$opt_gap||$opt_psi4||$opt_vasp||$opt_qe||$opt_python) {print "airss.pl: nbest not implemented\n";exit};

      my $spe=0.0;
      my $spe_best=1000000.0;
      my $k_best=0;

      for ($k = 1; $k <= $opt_nbest; $k++) {            
        system "buildcell < $opt_workdir/$uniquenew.cell > $opt_workdir/$unique.cell 2> /dev/null && cp $opt_workdir/$unique.cell $opt_workdir/$unique.cell.$k";
	if ( $opt_repose ) {
	  $spe=`(repose_relax 'repose -p $opt_pressure' $opt_mpinp $opt_workdir/$unique) | awk '/Enthalpy:/ {print \$2}'`;
	  system "cp $opt_workdir/$unique-out.cell $opt_workdir/$unique.cell.$k";
	  if ( $opt_castep ) {
	    system "cp $opt_workdir/$unique-out.cell $opt_workdir/$unique.cell";
	    $spe=`castep_relax -1 $executable $opt_sim $opt_symm $opt_workdir/$unique | awk '/Enthalpy:/ {print \$2}'`;
	  }
	} else {
	  $spe=`castep_relax -1 $executable $opt_sim $opt_symm $opt_workdir/$unique | awk '/Enthalpy:/ {print \$2}'`;
	}
        if ($spe < $spe_best) {$spe_best = $spe ; $k_best=$k};
      } 
      system "mv $opt_workdir/$unique.cell.$k_best $opt_workdir/$unique.cell ; rm -f $opt_workdir/$unique.cell.*";
    }

  } else {

    # Push a structure along a mode

    my $nions=`grep ions $uniquenew.phonon | awk '{printf \$4}'`;
    my $nbranches=`grep branches $uniquenew.phonon | awk '{printf \$4}'`;
  
    system "sed -n -e '1,/^%BLOCK positions/p' $uniquenew.cell > $unique.top";
    system "sed -n -e '/^%BLOCK [Pp][Oo][Ss]*/, /^%ENDBLOCK [Pp][Oo][Ss]*/p' $uniquenew.cell | sed '\$d' | sed '1d' > $unique.ions";
    system "sed -n -e '/^%ENDBLOCK [Pp][Oo][Ss]*/,\$p' $uniquenew.cell | sed -e '/^%BLOCK [Ee][Xx][Tt]*/, /^%ENDBLOCK [Ee][Xx][Tt]*/d' > $unique.bottom";
    
    my @latt=`sed -n -e '/^ Unit cell*/, /^ Frac*/p' $uniquenew.phonon | head -4 | tail -3`;
    
    my $nmode=int(($j+1)/2)+$opt_minmode-1;

    my $sign=(-1)**$j;
   
    my $ntail=$nions*($nbranches+1-$nmode);
    system "tail -$ntail $uniquenew.phonon | head -$nions | awk '{print \$3,\$5,\$7}' > $unique.mode";
    my @mode=`cat $unique.mode`;

	
    my $maxdist=0;
    foreach (@mode) {
      my $dist=0.0;
      my @vec = split(' ',$_);
      foreach (@latt) {
	my @latvec = split(' ',$_);
	$dist = $dist + ($latvec[0]*$vec[0]+$latvec[1]*$vec[1]+$latvec[2]*$vec[2])**2;
      }
      if (sqrt($dist) > $maxdist) {
	$maxdist=sqrt($dist);
      }
    }
   
    my $step=$sign*$opt_amp/$maxdist;

    system "cat $unique.top > $unique.cell";
    system "paste $unique.ions $unique.mode | awk '{print \$1,\$2+$step*\$5,\$3+$step*\$6,\$4+$step*\$7}' >> $unique.cell";
    system "cat $unique.bottom >> $unique.cell";
    system "rm -f $unique.top $unique.ions $unique.mode $unique.bottom";
   
  }
    
  # Add chosen external pressure
    
  open  CELLFILE, ">>$opt_workdir/$unique.cell" or die $!;
  print CELLFILE "%BLOCK EXTERNAL_PRESSURE\n";
  print CELLFILE  $opt_pressure," 0 0 \n";
  print CELLFILE  $opt_pressure," 0\n";
  print CELLFILE  $opt_pressure,"\n";
  print CELLFILE "%ENDBLOCK EXTERNAL_PRESSURE\n";
  close CELLFILE;
    
  # Perform the optimisation and store the results
    
  my @results='';

  if ( $opt_pp0 ) {
    @results=`pp0_relax pp0 $opt_workdir/$unique`;
  } elsif ( $opt_gosh )  {
    @results=`pp0_relax gosh $opt_workdir/$unique`;
  } elsif ( $opt_pp3 )  {
    @results=`pp3_relax pp3 $opt_workdir/$unique`;
  } elsif ( $opt_repose )  {
    @results=`(repose_relax repose $opt_mpinp $opt_workdir/$unique)`;
  } elsif ( $opt_build )  {
    @results="Volume: 0.00000001";
  } elsif ( $opt_gulp ) {
    if ( $opt_slab || $opt_cluster ) {
      @results=`gulp_relax gulp 1 $opt_pressure $opt_workdir/$unique`;
    } else {
      @results=`gulp_relax gulp 0 $opt_pressure $opt_workdir/$unique`;
    }
  } elsif ( $opt_lammps )  {
    @results=`lammps_relax lammps $opt_pressure $opt_workdir/$unique`;
  } elsif ( $opt_gap )  {
    @results=`gap_relax ase_relax_cell.py $opt_workdir/$unique`;
  } elsif ( $opt_psi4 )  {
    @results=`psi4_relax psi4 $opt_workdir/$unique`;
  } elsif ( $opt_vasp )  {
    @results=`vasp_relax $executable $opt_workdir/$unique`;
  } elsif ( $opt_qe ) {
    @results=`qe_relax $opt_steps $executable $opt_pressure $opt_sim $opt_symm $opt_workdir/$unique`;
  } elsif ( $opt_python ) {
    @results=`python_relax "python relax.py" $opt_pressure $opt_workdir/$unique`;
  } else {     
    @results=`castep_relax $opt_steps $executable $opt_sim $opt_symm $opt_workdir/$unique`;
  }

  # Perform an extra run of Castep if requested
  
  if ( $opt_castep && @results>0 ) {
    copy($opt_workdir."/".$unique."-out.cell",$opt_workdir."/".$unique.".cell");
    @results=`castep_relax $opt_steps $executable $opt_sim $opt_symm $opt_workdir/$unique`;
  }
 
  my $p=0;
  my $enthalpy=0;
  my $volume=0;
  foreach (@results) {
    my @tmp = split(' ',$_);
    if ($tmp[0] eq "Pressure:") {
      $p       =$tmp[1];
    }
    if ($tmp[0] eq "Enthalpy:") {
      $enthalpy=$tmp[1];
    }
    if ($tmp[0] eq "Volume:") {
      $volume  =$tmp[1];
    }
  }
     
  if ( ($enthalpy < $best_enthalpy) && ($volume > 0) ) {

    if ( ! $opt_track) {
      if (length($best_structure)>0) {
	system "rm -f $best_structure.* $opt_workdir/$best_structure.*";
      }
    }

    $best_structure = $unique;
	
    # If we are doing RASH

    if ($opt_num > 0) {
      $best_enthalpy = $enthalpy;
      $j=0;
    }

    # Calculate the DOS at the Fermi level

    if ($opt_dos) {
      system("(echo 'task : spectral' ; echo 'spectral_task : dos' ; grep -E -i -v 'task|write_bands|spin ' $opt_workdir/$unique.param ) > $opt_workdir/$unique-dos.param; cp $opt_workdir/$unique.cell $opt_workdir/$unique-dos.cell;eval $executable $opt_workdir/$unique-dos");
      system("(echo 'task : dos';echo 'broadening : linear'; echo 'compute_band_gap : true')>$opt_workdir/$unique-dos.odi;optados $opt_workdir/$unique-dos");
    }

    # Convert the structure to SHLX file, adding results, and computational parameters (do we need both versions?)
    
    if ($opt_nosymm == 0) {
      system("castep2res $opt_cluster $opt_workdir/$unique > $opt_workdir/$unique.res");
    } else {
      system("castep2res -1 $opt_workdir/$unique > $opt_workdir/$unique.res");
    }

    # Convert and store the results of the individual geometry optimisations
    
    if ($opt_harvest) {system("mkdir -p $opt_workdir/harvest && castepsplit2res < $opt_workdir/$unique.castep > $opt_workdir/harvest/$unique.res");}

    # Record the command line used
    
    system("sed -i 's/REM COMMAND_LINE/REM cmdline: airss.pl $commandline /g' $opt_workdir/$unique.res");
    
    # Identify and remove repeated structure, keeping track of the number of encounters (not needed if castep_relax is used)
    
    #if ($opt_sim > 0 && ($opt_pp3 || $opt_gulp) ) {
    if ($opt_sim > 0 ) {

      my $known ='';
      if ($opt_workdir eq ".") {
	$known = `find . -name "*.res" | xargs cat | cryan -c $opt_sim $unique 2> /dev/null | awk '{print \$1}' | head -1`;
      } else {
	# FIX THIS ...
	$known = '';
      }
      chomp($known);

      if ($known ne '') {
	my $knowntitl=`grep TITL $opt_workdir/$known.res | awk 'BEGIN {FS="n -"};{print \$1}'`;
	my $num_repeats=1+`grep TITL $opt_workdir/$known.res | awk 'BEGIN {FS="n -"};{print \$2}'`;
	chomp($knowntitl);

	open  RESFILE, "$opt_workdir/$known.res" or die $!;
	my @resdata = <RESFILE>;
	close RESFILE;
	open  RESFILE, ">$opt_workdir/$known.res" or die $!;
	shift(@resdata);
	unshift(@resdata,$knowntitl."n - ".$num_repeats."\n");
	print RESFILE @resdata;
	close RESFILE;

	if (! $opt_keep) {
	  system "rm -fr $opt_workdir/$unique*";
	} else {
	  system "mkdir -p $opt_workdir/trash ; mv $opt_workdir/$unique* $opt_workdir/trash";
	}
	next;
      }
      
    }

    # Calculate the vibrational modes

    if ( $opt_mode) {
      system "sed -e 's/geometryoptimization/thermodynamics/g' $opt_workdir/$unique.param > $opt_workdir/$unique.param.temp;
          echo 'phonon_method : finitedisplacement\nphonon_sum_rule : true\ncalculate_born_charges : false\nphonon_calc_lo_to_splitting : false' >> $opt_workdir/$unique.param.temp;
          echo 'PHONON_KPOINT_MP_GRID : 1 1 1' >> $opt_workdir/$unique.cell;
          mv $opt_workdir/$unique.param.temp $opt_workdir/$unique.param;
          eval $executable $opt_workdir/$unique;
          mv $opt_workdir/$unique.phonon $opt_workdir/$uniquenew.phonon";
    }

    # Make the current cell the new reference structure to shake
	
    updatepos ($opt_workdir."/".$uniquenew.".cell",$opt_workdir."/".$unique.".cell");
		
  } else {

    # Remove the output of the "failed" run
	
    if (! $opt_keep) {
      system "rm -fr $opt_workdir/$unique*";
      if ( $opt_qe ) {
         system "rm -fr $opt_workdir/CRASH";
         }
    } else {
      system "mkdir -p $opt_workdir/trash ; mv $opt_workdir/$unique* $opt_workdir/trash";
      if ( $opt_qe ) {
         system "if [ -f $opt_workdir/CRASH ]; then mv $opt_workdir/CRASH $opt_workdir/trash/$unique-CRASH; fi";
      }
    }
	
  }

  if ( $j < $opt_num ) {
    $j=$j+1;
    if ( $opt_amp < 0 ) {copy($opt_seed . ".cell",$opt_workdir."/".$uniquenew.".cell")};
  } else {
    $j=0;
    $best_enthalpy = $high_energy;
    $best_structure = "";
    copy($opt_seed . ".cell",$opt_workdir."/".$uniquenew.".cell");
  }

  if (! $opt_keep) {
    my @delete_files = <$opt_workdir/$unique*>;
    foreach (@delete_files) {if (($_ ne "$opt_workdir/$unique.res")&&($_ ne "$opt_workdir/$unique.conv")) {unlink($_);}}
    if ( $opt_vasp ) {
      my @delete_files = <$opt_workdir/$unique/*>;
      foreach (@delete_files) {unlink($_)};
      rmdir <$opt_workdir/$unique>;
    }
  }

  if ( $opt_best ) {
    my $best_formula='';
    if ( $opt_cluster ) {
      $best_formula=`cryan -r < $opt_workdir/$unique.res | awk '{print \$(NF-3),\$(NF-2)}'`;
    } else {
      $best_formula=`cryan -r < $opt_workdir/$unique.res | awk '{print \$(NF-2)}'`;
    }
    print $best_formula."\n";
    chomp($best_formula);    
    my $best_struct=`(find . -maxdepth 1 -name "*.res" ; echo $opt_workdir/$unique.res) | xargs cat | cryan -r | grep "$best_formula" | head -1 | awk '{print \$1}'`;
    chomp($best_struct);
    if ( $unique ne $best_struct ) {
      system "rm -fr $opt_workdir/$unique*";
    }
  }

  # Add this res file to the long list, and conv file if available
  
  if ( $opt_pack ) {
    system "cat $opt_workdir/$unique.res >> $opt_workdir/$uniquenew.res && rm -f $opt_workdir/$unique.res";
    if ( -e "$opt_workdir/$unique.conv" ) {
      system "cat $opt_workdir/$unique.conv >> $opt_workdir/$uniquenew.conv && rm -f $opt_workdir/$unique.conv"
    }
  } 
  
  # Use rsync to copy files from the work directory

  if ($opt_workdir ne ".") {
    system "rsync -t $opt_workdir/$uniquehead* .";
  }
  
}

sub updatepos {

  my $oldcell=shift;
  my $newcell=shift;

  system "sed -n -e '/^%BLOCK [Ll][Aa][Tt]*/, /^%ENDBLOCK [Pp][Oo][Ss]*/p' $newcell > $newcell.temp;
    sed -e '/ANG/d' $newcell.temp > $oldcell.temp;
    sed -e '/^%BLOCK [Ll][Aa][Tt]*/, /^%ENDBLOCK [Ll][Aa][Tt]*/d' $oldcell | sed -e '/^%BLOCK [Pp][Oo][Ss]*/, /^%ENDBLOCK [Pp][Oo][Ss]*/d' >> $oldcell.temp;
    sed -e 's/latTICE/LATTICE/g' $oldcell.temp | grep -v '\#VARVOL=' | grep -v '\#TARGVOL=' | grep -v '\#POSAMP=' | grep -v '\#CELLAMP=' | grep -v '\#NATOM=' | grep -v '\#SPECIES=' | grep -v '\#FOCUS=' | grep -v '\#SYMM' | grep -v '\#NFORM=' | grep -v '\#SUPER' > $oldcell;
    echo '\#POSAMP=$opt_amp' >> $oldcell;
    echo '\#CELLAMP=$opt_camp' >> $oldcell;
    rm $newcell.temp $oldcell.temp;";
}