/*
Copyright (C) 2003 by Sean David Fleming

sean@ivec.org

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

The GNU GPL can also be found at http://www.gnu.org
*/

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>

#include "gdis.h"
#include "coords.h"
#include "edit.h"
#include "file.h"
#include "gulp_keyword.h"
#include "parse.h"
#include "scan.h"
#include "space.h"
#include "matrix.h"
#include "model.h"
#include "numeric.h"
#include "interface.h"

#define DEBUG_MORE 0
#define MAX_KEYS 15

/* main structures */
extern struct sysenv_pak sysenv;
extern struct elem_pak elements[];


/*********************************/
/* GULP structure initialization */
/*********************************/
/* TODO - put everything GULP related from model.c in here */
void gulp_init(gpointer data)
{
struct gulp_pak *gulp = data;

g_assert(gulp != NULL);

gulp->solvation_model = GULP_SOLVATION_NONE;
gulp->cosmo_sas = TRUE;
gulp->cosmo_shape = 0;
gulp->cosmo_shape_index[0] = 5;
gulp->cosmo_shape_index[1] = 0;
gulp->cosmo_points = 974;
gulp->cosmo_segments = 110;
gulp->cosmo_smoothing = 0.0;
gulp->cosmo_solvent_epsilon = 1.0;
gulp->cosmo_solvent_radius = 1.0;
gulp->cosmo_solvent_delta = 0.1;
gulp->cosmo_solvent_rmax = 10.0;
}

/***********************************************************/
/* use number of points and shape to compute shape indices */
/***********************************************************/
#define COSMO_SHAPE_INDEX_MAX 10
gint gulp_shape_indices_compute(struct gulp_pak *gulp)
{
gint i, j, n;
gdouble a, b, c;

g_assert(gulp != NULL);

switch (gulp->cosmo_shape)
  {
  case 0:
    c = 4.0;
    break;
  default:
    c = 10.0;
    break;
  }

a = gulp->cosmo_points;
a = log((a - 2.0) / c) / log(3.0);
b = log(4.0) / log(3.0);

/* simple for loop search for suitable index values */
for (j=0 ; j<COSMO_SHAPE_INDEX_MAX ; j++)
  {
  i = nearest_int(a - j*b);

  n = nearest_int (c * pow(3.0, (gdouble) i) * pow(4.0, (gdouble) j) + 2.0);

/*
printf("[%d,%d] : %d / %d\n", i, j, n, gulp->cosmo_points);
*/

  if (!(n - gulp->cosmo_points))
    {
    gulp->cosmo_shape_index[0] = i;
    gulp->cosmo_shape_index[1] = j;

    return(TRUE);
    }
  }
return(FALSE);
}

/*****************************/
/* core file print primitive */
/*****************************/
void fprint_core(FILE *fp, struct core_pak *core, struct model_pak *model)
{
gint i;
gdouble vec[3];

/* get cartesian coords */
ARR3SET(vec, core->x);
vecmat(model->latmat, vec); 

/* set fractional part */
for (i=0 ; i<model->periodic ; i++)
  if (model->fractional)
    vec[i] = core->x[i];

/* only print specific atom charges if we read them in */
/* in that fashion (ie lookup_charge is false) */
if (core->breathe)
  fprintf(fp,"%-4s bcore %11.6f %11.6f %11.6f", core->atom_label,
                                                vec[0], vec[1], vec[2]);
else
  fprintf(fp,"%-4s core  %11.6f %11.6f %11.6f", core->atom_label,
                                                vec[0], vec[1], vec[2]);

/* NB: 7 dp for charges - 8 occasionally introduced small errors (cf GULP) */
if (model->gulp.print_charge)
  if (!core->lookup_charge)
    fprintf(fp,"  %11.7lf", core->charge);

/* 7dp for occupancy - to cover cases like 1/3 */
if (core->has_sof)
  fprintf(fp,"  %9.7f", core->sof);

if (core->flags)
  fprintf(fp," %s", core->flags);

if (model->periodic == 2)
  if (core->growth)
    fprintf(fp, " %%");

if (core->translate)
  fprintf(fp, " T");

fprintf(fp,"\n");
}

/******************************/
/* shell file print primitive */
/******************************/
void fprint_shell(FILE *fp, struct shel_pak *shell, struct model_pak *model)
{
gint i;
gdouble vec[3];

/* get cartesian coordinates */
ARR3SET(vec, shell->x);
vecmat(model->latmat, vec); 

/* set fractional part */
for (i=0 ; i<model->periodic ; i++)
  if (model->fractional)
    vec[i] = shell->x[i];

if (shell->breathe)
  fprintf(fp,"%-4s bshel %11.6f %11.6f %11.6f", shell->shell_label,
                                                vec[0], vec[1], vec[2]);
else
  fprintf(fp,"%-4s shel  %11.6f %11.6f %11.6f", shell->shell_label,
                                                vec[0], vec[1], vec[2]);

/* only print specific atom charges if we read them in */
if (model->gulp.print_charge)
  if (!shell->lookup_charge)
    fprintf(fp,"  %11.7lf", shell->charge);

/* 7dp for occupancy - to cover cases like 1/3 */
if (shell->has_sof)
  fprintf(fp,"  %9.7f", shell->sof);

if (shell->flags)
  fprintf(fp," %s", shell->flags);

if (model->periodic == 2)
  if ((shell->core)->growth)
    fprintf(fp, " %%");

if (shell->translate)
  fprintf(fp, " T");

fprintf(fp,"\n");
}

/*********************************************/
/* enforce valid input/output GULP filenames */
/*********************************************/
void gulp_files_init(struct model_pak *model)
{
gchar *basename;

g_assert(model != NULL);
g_assert(model->basename != NULL);

/* remove spaces from basename */
basename = g_strdup(model->basename);
parse_space_replace(basename, '_');

/* must have input and ouptut filenames */
if (!model->gulp.temp_file)
  model->gulp.temp_file = g_strdup_printf("%s.gin", model->basename);
if (!model->gulp.out_file)
  model->gulp.out_file = g_strdup_printf("%s.got", model->basename);
/* default dump file - gulp opti ASSUMES a dump file - see proc_gulp_task() */
if (!model->gulp.dump_file)
  model->gulp.dump_file = g_strdup_printf("%s.res", model->basename);
}

/*********************************/
/* cleanup for GULP file control */
/*********************************/
void gulp_files_free(struct model_pak *model)
{
g_assert(model != NULL);

g_free(model->gulp.temp_file);
g_free(model->gulp.dump_file);
g_free(model->gulp.out_file);
}

/**********************************/
/* free memory used for gulp data */
/**********************************/
#define DEBUG_FREE_GULP_DATA 0
void gulp_data_free(struct model_pak *model)
{
g_assert(model != NULL);

/* NEW - GULP files with associated trajectory only */
if (model->gulp.orig_cores)
  free_slist(model->gulp.orig_cores);
if (model->gulp.orig_shells)
  free_slist(model->gulp.orig_shells);

g_free(model->gulp.potentials);
g_free(model->gulp.elements);
g_free(model->gulp.species);
g_free(model->gulp.kpoints);
model->gulp.potentials=NULL;
model->gulp.elements=NULL;
model->gulp.species=NULL;
model->gulp.kpoints=NULL;

if (model->gulp.epot_vecs)
  free_slist(model->gulp.epot_vecs);
if (model->gulp.epot_vals)
  free_slist(model->gulp.epot_vals);
model->gulp.epot_vecs=NULL;
model->gulp.epot_vals=NULL;

g_free(model->gulp.libfile);
g_free(model->gulp.extra);
g_free(model->gulp.eatt_units);
g_free(model->gulp.esurf_units);
g_free(model->gulp.pressure);
g_free(model->gulp.temperature);
model->gulp.libfile=NULL;
model->gulp.extra=NULL;
model->gulp.eatt_units=NULL;
model->gulp.esurf_units=NULL;
model->gulp.pressure=NULL;
model->gulp.temperature=NULL;
}

/***********************************************/
/* copy unrecognized gulp keywords and options */
/***********************************************/
void gulp_extra_copy(struct model_pak *src, struct model_pak *dest)
{
if (src->gulp.extra)
  {
  g_free(dest->gulp.extra);
  dest->gulp.extra = g_strdup(src->gulp.extra);
  }
if (src->gulp.extra_keywords)
  {
  g_free(dest->gulp.extra_keywords);
  dest->gulp.extra_keywords = g_strdup(src->gulp.extra_keywords);
  }
}

/**********************/
/* transfer gulp data */
/**********************/
void gulp_data_copy(struct model_pak *src, struct model_pak *dest)
{
/* free any old data in the destination model */
gulp_data_free(dest);

/* copy keywords */
dest->gulp.no_exec = src->gulp.no_exec;
dest->gulp.run = src->gulp.run;
dest->gulp.method = src->gulp.method;
dest->gulp.optimiser = src->gulp.optimiser;
dest->gulp.optimiser2 = src->gulp.optimiser2;
dest->gulp.switch_type = src->gulp.switch_type;
dest->gulp.switch_value = src->gulp.switch_value;
dest->gulp.free = src->gulp.free;
dest->gulp.qeq = src->gulp.qeq;
dest->gulp.zsisa = src->gulp.zsisa;
dest->gulp.compare = src->gulp.compare;
dest->gulp.nosym = src->gulp.nosym;
dest->gulp.fix = src->gulp.fix;
dest->gulp.cycles = src->gulp.cycles;
dest->gulp.ensemble = src->gulp.ensemble;
dest->gulp.coulomb = src->gulp.coulomb;
dest->gulp.maxcyc = src->gulp.maxcyc;

/* FIXME AR change - copy over cosmo stuff - may cause problems with surfaces (but thats what I want) */
/* Note to SEAN would be a lot easier if were a separate struct as I could then replace with a single memcopy */
dest->gulp.solvation_model = src->gulp.solvation_model;
dest->gulp.cosmo_sas = src->gulp.cosmo_sas; 
dest->gulp.cosmo_shape = src->gulp.cosmo_shape; 
dest->gulp.cosmo_shape_index[0] = src->gulp.cosmo_shape_index[0]; 
dest->gulp.cosmo_shape_index[1] = src->gulp.cosmo_shape_index[1];
dest->gulp.cosmo_points = src->gulp.cosmo_points;
dest->gulp.cosmo_segments = src->gulp.cosmo_segments;
dest->gulp.cosmo_smoothing = src->gulp.cosmo_smoothing; 
dest->gulp.cosmo_solvent_epsilon = src->gulp.cosmo_solvent_epsilon;
dest->gulp.cosmo_solvent_radius = src->gulp.cosmo_solvent_radius;
dest->gulp.cosmo_solvent_delta = src->gulp.cosmo_solvent_delta;
dest->gulp.cosmo_solvent_rmax = src->gulp.cosmo_solvent_rmax;

/* FIXME - may cause problems when creating surfaces from unit cells */
dest->gulp.output_extra_keywords = src->gulp.output_extra_keywords;
dest->gulp.output_extra = src->gulp.output_extra;
gulp_extra_copy(src, dest);

/* copy surface related data */
if (dest->periodic == 2)
  dest->gulp.sbulkenergy = src->gulp.sbulkenergy;
if (dest->periodic == 2 && src->periodic == 2)
  dest->surface.dspacing = src->surface.dspacing;

/* duplicate gulp character data */
dest->gulp.libfile = g_strdup(src->gulp.libfile);
dest->gulp.potentials = g_strdup(src->gulp.potentials);
dest->gulp.elements = g_strdup(src->gulp.elements);
dest->gulp.species = g_strdup(src->gulp.species);
dest->gulp.pressure = g_strdup(src->gulp.pressure);
dest->gulp.temperature = g_strdup(src->gulp.temperature);
}

/**********************************/
/* debugging - dump a phonon mode */
/**********************************/
/* NB: modes are from 1 - num_phonons */
void dump_phonon(gint atom, gint mode, struct model_pak *model)
{
gchar *freq;
gdouble v[3];
GSList *xlist, *ylist, *zlist;
struct core_pak *core;

freq = g_slist_nth_data(model->phonons, mode-1);
core = g_slist_nth_data(model->cores, atom-1);

xlist = core->vibx_list; 
ylist = core->viby_list; 
zlist = core->vibz_list; 

v[0] = *((gdouble *) g_slist_nth_data(xlist, mode-1));
v[1] = *((gdouble *) g_slist_nth_data(ylist, mode-1));
v[2] = *((gdouble *) g_slist_nth_data(zlist, mode-1));

printf("[%s] : ", freq);
printf("[%s]", core->atom_label);
P3VEC(" : ", v);
}

/************************************************/
/* debugging - dump the vibrational frequencies */
/************************************************/
void dump_phonons(struct model_pak *model)
{
GSList *list;

for (list=model->phonons ; list ; list=g_slist_next(list))
  {
  printf("[%s]\n", (gchar *) list->data);
  }
}

/**********************************/
/* print appropriate region label */
/**********************************/
void write_region_header(FILE *fp, gint region, struct model_pak *model)
{
/* check periodicity */
switch (model->periodic)
  {
  case 1:
    if (model->fractional)
      fprintf(fp,"pfractional region %d", region+1);
    else
      fprintf(fp,"cartesian region %d", region+1);
    break;

  case 2:
    if (model->fractional)
      fprintf(fp,"sfractional region %d", region+1);
    else
      fprintf(fp,"cartesian region %d", region+1);
    break;

  case 3:
    if (model->fractional)
      fprintf(fp,"fractional");
    else
      fprintf(fp,"cartesian");
    break;

  default:
    fprintf(fp,"cartesian");
  }

/* NB: rigid xyz allows all translation as a rigid body */
/* rigid z allows vertical movement only -> mandatory for grid based dock */
if (model->gulp.rigid && region == 2)
  {
  fprintf(fp," rigid ");
  if (model->gulp.rigid_move)
    fprintf(fp,"%s ", model->gulp.rigid_move);
  }

fprintf(fp,"\n");
}

/*********************/
/* write a GULP file */
/*********************/
#define DEBUG_WRITE_GULP 0
gint write_gulp(gchar *filename, struct model_pak *data)
{
gint libflag, flag, region;
GString *line;
GSList *list;
struct core_pak *core;
struct shel_pak *shell;
struct vec_pak *vec;
FILE *fp;

g_return_val_if_fail(data != NULL, 1);

/* is this really necessary??? */
unlink(filename);

/* is a valid potential library specified? */
libflag = 0;
if (data->gulp.libfile)
  {
  if (strlen(data->gulp.libfile))
    {
/* test for library existance */
    if (!access(data->gulp.libfile, R_OK))
      libflag++;
    else
      printf("Warning: library %s not found.\n", data->gulp.libfile);
    }
  }

/* off we go */
fp = fopen(filename, "wt");
if (!fp)
  {
  printf("Failed to open %s.\n", filename);
  return(2);
  }

/* run type */
line = g_string_new(NULL);
switch(data->gulp.run)
  {
  case E_SINGLE:
    g_string_assign(line,"single");
    break;
  case E_OPTIMIZE:
    g_string_assign(line,"opti");
    break;
  case MD:
    g_string_assign(line,"md");
    break;
/* no run -> single pt (valid default?) */
  default:
    g_string_assign(line,"single");
    break;
  }

switch(data->gulp.method)
  {
  case CONP:
    g_string_append(line," conp");
    break;
  case CONV:
    g_string_append(line," conv");
    break;
  default:
    break;
  }

switch(data->gulp.optimiser)
  {
  case BFGS_OPT:
    g_string_append(line," bfgs");
    break;
  case CONJ_OPT:
    g_string_append(line," conjugate");
    break;
  case RFO_OPT:
    g_string_append(line," rfo");
    break;
  }

if (data->gulp.unit_hessian)
    g_string_append(line," unit");

switch(data->gulp.coulomb)
  {
  case MOLE:
    g_string_append(line," mole");
    break;
  case MOLMEC:
    g_string_append(line," molmec");
    break;
  case MOLQ:
    g_string_append(line," molq");
    break;
  default:
    break;
  }

/* cosmo calc */
/* TODO - change variable name to solvation model */
switch (data->gulp.solvation_model)
  {
  case GULP_SOLVATION_COSMIC:
    g_string_append(line," cosmic");
    break;
  case GULP_SOLVATION_COSMO:
    g_string_append(line," cosmo");
    break;
  }

/* phonon calc */
if (data->gulp.phonon)
  {
  g_string_append(line," phonon");
  if (data->gulp.eigen)
    g_string_append(line," eigen");
  }

/* NEW - this makes it easy to use GULP to compute properties */
if (data->gulp.prop)
  g_string_append(line," eigen");

/* additional keywords */
if (data->gulp.qeq)
  g_string_append(line," qeq");
if (data->gulp.free)
  g_string_append(line," free");
if (data->gulp.zsisa)
  g_string_append(line," zsisa");
if (data->gulp.compare)
  g_string_append(line," comp");
if (data->gulp.fix)
  g_string_append(line," fix");
/* if nosym specified, force output of non-prim cell for consistency */
if (data->gulp.nosym)
  g_string_append(line," nosym full");
/* electrostatic potential calcs */
if (g_slist_length(data->gulp.epot_vecs))
  g_string_append(line," epot");

/* print out the keywords line */
if (data->gulp.output_extra_keywords && data->gulp.extra_keywords)
  fprintf(fp,"%s %s\n\n", line->str, data->gulp.extra_keywords);
else
  fprintf(fp,"%s\n\n",line->str);

fprintf(fp, "# ");
gdis_blurb(fp);
fprintf(fp, "#\n\n");

if (data->periodic == 2)
  {
/* NB: this surface data is only meaningful if generate from within GDIS */
  if (VEC3MAGSQ(data->surface.depth_vec) > FRACTION_TOLERANCE)
    {
    fprintf(fp, "#       miller: %d %d %d\n",
                     (gint) data->surface.miller[0],
                     (gint) data->surface.miller[1],
                     (gint) data->surface.miller[2]);
    fprintf(fp, "#        shift: %lf\n", data->surface.shift);
    fprintf(fp, "# region sizes: %d %d\n\n",
                    (gint) data->surface.region[0],
                    (gint) data->surface.region[1]);
    }
  }

/* optional data */
if (data->gulp.maxcyc > 0)
  fprintf(fp,"maxcyc %d\n\n", (gint) data->gulp.maxcyc);

if (data->gulp.optimiser2 != SWITCH_OFF)
  {
  fprintf(fp,"switch ");
  switch(data->gulp.optimiser2)
    {
    case CONJ_OPT:
      fprintf(fp,"conj ");
      break;
    case RFO_OPT:
      fprintf(fp,"rfo ");
      break;
    case BFGS_OPT:
    default:
      fprintf(fp,"bfgs ");
      break;
    }

  switch(data->gulp.switch_type)
    {
    case CYCLE:
      fprintf(fp,"cycle %d\n\n", (gint) data->gulp.switch_value);
      break;
    case GNORM:
    default:
      fprintf(fp,"gnorm %lf\n\n", data->gulp.switch_value);
      break;
    }
  }

/* only write if we have a valid supplied library */
if (libflag)
  fprintf(fp,"library %s\n\n", data->gulp.libfile);

fprintf(fp,"name %s\n\n", data->basename);

if (data->gulp.dump_file)
  if (strlen(data->gulp.dump_file))
    fprintf(fp,"\ndump %s\n\n", data->gulp.dump_file);

if (data->gulp.temperature)
  if (strlen(data->gulp.temperature))
    fprintf(fp,"temperature %s\n\n", data->gulp.temperature);

if (data->gulp.pressure)
  if (strlen(data->gulp.pressure))
    fprintf(fp,"pressure %s\n\n", data->gulp.pressure);

/* dynamics stuff */
if (data->gulp.run == MD)
  {
  switch (data->gulp.ensemble)
    {
    case NVE:
      fprintf(fp,"ensemble nve\n\n");
      break;
    case NVT:
      fprintf(fp,"ensemble nvt\n\n");
      break;
    case NPT:
      fprintf(fp,"ensemble npt\n\n");
      break;
    }

  fprintf(fp, "timestep %s\n", data->gulp.timestep);
  fprintf(fp, "equilibration %s\n", data->gulp.equilibration);
  fprintf(fp, "production %s\n", data->gulp.production);
  fprintf(fp, "sample %s\n", data->gulp.sample);
  fprintf(fp, "write %s\n", data->gulp.write);
  fprintf(fp, "\n");

  if (data->gulp.trj_file)
    fprintf(fp, "output trajectory %s\n", data->gulp.trj_file);
  }

if (data->gulp.mov_file)
  {
  fprintf(fp, "output movie arc %s\n", data->gulp.mov_file);
  if (data->gulp.run == MD)
    fprintf(fp, "mdarchive %s\n", data->gulp.mov_file);
  }


/* NEW - cosmo stuff */
if (data->gulp.solvation_model != GULP_SOLVATION_NONE)
  {
  switch (data->gulp.cosmo_shape)
    {
    case 0:
      fprintf(fp, "cosmoshape octahedron\n");
      break;
    default:
      fprintf(fp, "cosmoshape dodecahedron\n");
    }
  fprintf(fp, "pointsperatom %d\n", data->gulp.cosmo_points);
  fprintf(fp, "segmentsperatom %d\n", data->gulp.cosmo_segments);

  fprintf(fp, "solventepsilon %f\n", data->gulp.cosmo_solvent_epsilon);
  fprintf(fp, "solventradius %f %f\n", data->gulp.cosmo_solvent_radius, data->gulp.cosmo_solvent_delta);
  fprintf(fp, "solventrmax %f\n", data->gulp.cosmo_solvent_rmax);
  fprintf(fp, "rangeforsmooth %f\n", data->gulp.cosmo_smoothing);

  fprintf(fp, "\n");

  if (data->gulp.cosmo_sas)
    fprintf(fp, "output sas %s.sas\n", data->basename);
  }


#if DEBUG_WRITE_GULP
printf("  periodic: %d\n", data->periodic);
printf("fractional: %d\n", data->fractional);
#endif

switch(data->periodic)
  {
  case 3:
    if (data->construct_pbc)
      {
/* NB: gulp matrices are transposed wrt gdis */
      fprintf(fp,"vectors\n");
      fprintf(fp,"%15.10f %15.10f %15.10f\n",
                  data->latmat[0], data->latmat[3], data->latmat[6]);
      fprintf(fp,"%15.10f %15.10f %15.10f\n",
                  data->latmat[1], data->latmat[4], data->latmat[7]);
      fprintf(fp,"%15.10f %15.10f %15.10f\n",
                  data->latmat[2], data->latmat[5], data->latmat[8]);
      }
    else
      {
      fprintf(fp,"cell\n");
      fprintf(fp,"%lf %lf %lf  %lf %lf %lf\n",data->pbc[0],data->pbc[1],data->pbc[2],
                data->pbc[3]*180/PI, data->pbc[4]*180/PI, data->pbc[5]*180/PI);
      }
    break;

  case 2:
    if (data->construct_pbc)
      {
      fprintf(fp,"svectors\n");
      fprintf(fp,"%lf  %lf\n%lf  %lf\n",data->latmat[0], data->latmat[3],
                                    data->latmat[1], data->latmat[4]);
      }
    else
      {
      fprintf(fp,"scell\n");
      fprintf(fp,"%lf %lf %lf\n",data->pbc[0], data->pbc[1],
                              180.0*data->pbc[5]/PI);
      }
    break;

  case 1:
    if (data->construct_pbc)
      {
      fprintf(fp,"pvector\n");
      fprintf(fp,"%lf\n",data->latmat[0]);
      }
    else
      {
      fprintf(fp,"pcell\n");
      fprintf(fp,"%lf\n",data->pbc[0]);
      }
    break;
  }

/* coord section write */
flag = 0;

/* CURRENT - no more than 3 regions allowed */
for (region=0 ; region<3 ; region++)
  {
/* cores */
  for (list=data->cores ; list ; list=g_slist_next(list))
    {
    core = list->data;

    if (core->status & DELETED)
      continue;
    if (!core->primary)
      continue;
    if (core->region != region)
       continue;

/* write suitable header if it hasn't yet been written */
    if (!flag)
      {
      write_region_header(fp, region, data);
      flag++;
      }
    fprint_core(fp, core, data);
    }

/* shells */
  for (list=data->shels ; list ; list=g_slist_next(list))
    {
    shell = list->data;

    if (shell->status & DELETED)
      continue;
    if (!shell->primary)
      continue;
    if (shell->region != region)
      continue;

    fprint_shell(fp, shell, data);
    }
  flag = 0;
  }

/* post coord data */
switch(data->periodic)
  {
  case 2:
/* write dhkl for Eatt calcs */
    if (data->surface.dspacing > 0.1 && !data->gulp.no_eatt)
      fprintf(fp, "\ndhkl %lf\n",data->surface.dspacing);

/* write surface bulk energy for Esurf calcs */
    if (!data->gulp.no_esurf)
      fprintf(fp, "\nsbulkenergy %lf\n",data->gulp.sbulkenergy);
    break;

/* SG info */
  case 3:
    if (data->sginfo.spacename)
      {
      fprintf(fp,"\nspace\n%s\n",g_strstrip(data->sginfo.spacename));
      if (data->sginfo.cellchoice)
        fprintf(fp,"\norigin %d\n",data->sginfo.cellchoice);
      }
    else
      {
      fprintf(fp,"\nspace\nP 1\n");
      }
    break;
  }
fprintf(fp,"\n");

/* output species data */
if (data->gulp.species)
  {
  strip_extra(data->gulp.species);
  if (strlen(data->gulp.species))
    {
    fprintf(fp,"species\n");
    fprintf(fp,"%s\n", data->gulp.species);
    fprintf(fp,"end\n\n");
    }
  }

/* output cova data */
/* FIXME - what if other element properties have been modified (within GDIS) */
if (data->gulp.elements)
  {
  strip_extra(data->gulp.elements);
  if (strlen(data->gulp.elements))
    {
    fprintf(fp,"elements\n");
    fprintf(fp,"%s\n", data->gulp.elements);
    fprintf(fp,"end\n\n");
    }
  }

if (data->gulp.potentials)
  if (strlen(data->gulp.potentials))
  fprintf(fp,"%s\n", data->gulp.potentials);

if (data->gulp.epot_vecs)
  {
  fprintf(fp, "potsites cart\n");

  for (list=data->gulp.epot_vecs ; list ; list=g_slist_next(list))
    {
    vec = list->data;

    fprintf(fp, "%f  %f  %f\n", vec->rx[0], vec->rx[1], vec->rx[2]);
    }
  fprintf(fp,"\n");
  }

/* specified kpoints? */
if (data->periodic && data->gulp.kpoints)
  {
  strip_extra(data->gulp.kpoints);
  if (strlen(data->gulp.kpoints))
    {
    fprintf(fp, "kpoints %d\n", 1+char_count(data->gulp.kpoints, '\n'));
    fprintf(fp, "%s\n\n", data->gulp.kpoints);
    }
  }

/* unrecognized stuff */
if (data->gulp.output_extra)
  {
  if (data->gulp.extra)
    {
    fprintf(fp, "\n#--- Additional options/data\n\n");
    fprintf(fp, "%s\n\n", data->gulp.extra);
    }
  }
fprintf(fp, "print 1\n");

/* done */
fclose(fp);
g_string_free(line, TRUE);
return(0);
}

/*****************/
/* debugging aid */
/*****************/
void print_gulp_flags(gint *list, gint size)
{
gint i;

printf("Flags:");
for (i=0 ; i<size ; i++)
  {
/* active? */
  if (*(list+i))
    {
    switch(i)
      {
      case NAME:
        printf(" name");
        break;
      case SURFACE_CELL:
        printf(" scell");
        break;
      case CELL:
        printf(" cell");
        break;
      case FRAC:
        printf(" frac");
        break;
      case SFRAC:
        printf(" sfrac");
        break;
      case PFRAC:
        printf(" pfrac");
        break;
      case CART:
        printf(" cart");
        break;
      case DUMP_FILE:
        printf(" dump");
        break;
      case SURFACE_VECTORS:
        printf(" svec");
        break;
      case LATTICE_VECTORS:
        printf(" lvec");
        break;
      }
    }
  }
printf("\n");
}

/***********************/
/* gulp keyword parser */
/***********************/
gint gulp_is_keyword(gchar *text)
{
gint i, n;

if (!text)
  return(-1);

i=0;
while (gulp_keyword[i].label)
  {
  n = strlen(gulp_keyword[i].label);
  if (g_strncasecmp(text, gulp_keyword[i].label, n) == 0)
    return(gulp_keyword[i].code);
  i++;
  }
return(-1);
}

/***********************/
/* gulp options parser */
/***********************/
gint gulp_is_option(gchar *text)
{
gint i, n;

if (!text)
  return(-1);

i=0;
while (gulp_option[i].label)
  {
  n = strlen(gulp_option[i].label);
  if (g_strncasecmp(text, gulp_option[i].label, n) == 0)
    return(gulp_option[i].code);
  i++;
  }
return(-1);
}

enum
{
GULP_NEW_NAME, GULP_NEW_LATTICE, GULP_NEW_COORDS, GULP_NEW_FLAGS
};

/************************************************/
/* determines if we should allocate a new model */
/************************************************/
struct model_pak *gulp_change_model(gint trigger, struct model_pak *model)
{
gint i, flag_new=FALSE;
static gint flag_level=1, flag_table[GULP_NEW_FLAGS];
struct model_pak *new;

if (!model)
  {
/* init flag table */
  flag_level = 1;
  for (i=GULP_NEW_FLAGS ; i-- ; )
    flag_table[i] = 0;
  return(NULL);
  }

/* handle possible new model triggers (NB: could occur in any order) */
switch (trigger)
  {
  case GULP_NEW_COORDS:
  case GULP_NEW_NAME:
  case GULP_NEW_LATTICE:
    flag_table[trigger]++;
/* only new flag triggers above the current level trigger a new model */
    if (flag_table[trigger] > flag_level)
      {
      flag_new = TRUE;
      flag_level = flag_table[trigger];
      }
    break;
  }

/* once the new coords flag is tripped - any subsequent new flag */
/* triggers should generate a new model. */
/* NB: name/cell options must come before the coords (GULP incompatibility?) */
if (trigger == GULP_NEW_COORDS)
  {
  for (i=GULP_NEW_FLAGS ; i-- ; )
    flag_table[i] = flag_level;
  }

/*
printf("%d : %d : %d\n", flag_table[0], flag_table[1], flag_table[2]); 
*/

/* return new model if flagged, else the old model */
if (flag_new)
  {
/*
printf(" + new model triggered!\n");
*/
  new = model_new();
  new->id = -1;
  return(new);
  }
return(model);
}

/*********************************/
/* yet another gulp read rewrite */
/*********************************/
#define DEBUG_READ_GULP 0
gint read_gulp(gchar *filename, struct model_pak *model)
{
gint i, j, nc, ns;
gint keywords_found, keywords_expected, num_tokens, remaining;
gint code, context, run, method, optimiser, optimiser2, unit, fix, switch_type;
gint region, breathe, coulomb, maxcyc, qeq, free, zsisa, compare, nosym, phonon, eigen;
gint type, gen_flag, temp_code, noflags, growth, translate;
gchar *tmp, *line, **buff;
gdouble switch_value;
gdouble vec1[3], vec2[3], vec3[3];
gpointer scan;
GSList *item, *clist, *slist;
GString *key_buff, *elem_buff, *ex_buff, *ff_buff, *species_buff, *kpoints_buff;
struct gulp_pak gulp;
struct elem_pak elem_data;
struct core_pak *core;
struct shel_pak *shell;

scan = scan_new(filename);
if (!scan)
  return(1);

/* init */
model->id = -1;
keywords_found = keywords_expected = 0;
/* deprec - use gulp_pak structure/gulp_init() instead */
code = context = run = method = optimiser = optimiser2 = switch_type = switch_value = -1;
qeq = free = zsisa = compare = nosym = phonon = eigen = unit = fix = FALSE;

/* NEW */
gulp_init(&gulp);

/* NEW - to fix GULP's out of control coord lines - assume noflags */
noflags = TRUE;
maxcyc = 0;
region = REGION1A;
coulomb = NOBUILD;
gulp_change_model(0, NULL);
key_buff = g_string_new(NULL);
elem_buff = g_string_new(NULL);
ex_buff = g_string_new(NULL);
ff_buff = g_string_new(NULL);
species_buff = g_string_new(NULL);
kpoints_buff = g_string_new(NULL);

/* main parse loop */
buff = scan_get_tokens(scan, &num_tokens);
while (!scan_complete(scan))
  {

g_assert(buff != NULL);

  if (!keywords_found)
    {
/* keyword only search */
    for (i=num_tokens ; i-- ; )
      {
      code = gulp_is_keyword(*(buff+i));
      switch (code)
        {
        case E_SINGLE:
        case E_OPTIMIZE:
        case MD:
          run = code;
          keywords_found++;
          keywords_expected = num_tokens;
          break;
        case RFO_OPT:
        case BFGS_OPT:
        case CONJ_OPT:
          optimiser = code;
          keywords_found++;
          keywords_expected = num_tokens;
          break;
        case UNIT_HESSIAN:
          unit = TRUE;
          keywords_found++;
          keywords_expected = num_tokens;
          break;
        case FIX:
          fix = TRUE;
          keywords_found++;
          break;
        case MOLE:
        case MOLMEC:
        case MOLQ:
          coulomb = code;
          keywords_found++;
          keywords_expected = num_tokens;
          break;
        case QEQ:
          qeq = TRUE;
          keywords_found++;
          keywords_expected = num_tokens;
          break;
        case FREE_ENERGY:
          free = TRUE;
          keywords_found++;
          keywords_expected = num_tokens;
          break;
        case ZSISA:
          zsisa = TRUE;
          keywords_found++;
          keywords_expected = num_tokens;
          break;
        case COMPARE:
          compare = TRUE;
          keywords_found++;
          keywords_expected = num_tokens;
          break;
        case NOSYM:
          nosym = TRUE;
          keywords_found++;
          keywords_expected = num_tokens;
          break;
        case CONP:
        case CONV:
          method = code;
          keywords_found++;
          keywords_expected = num_tokens;
          break;

/*
        case NOFLAGS:
          noflags = TRUE;
          break;
*/


        case GULP_SOLVATION_COSMIC:
        case GULP_SOLVATION_COSMO:
          gulp.solvation_model = code;
          keywords_found++;
          keywords_expected = num_tokens;
          break;


        case PHONON:
          phonon = TRUE;
          keywords_found++;
          keywords_expected = num_tokens;
          break;
        case EIGEN:
          eigen = TRUE;
          keywords_found++;
          keywords_expected = num_tokens;
          break;

        default:
          g_string_sprintfa(key_buff, "%s ", *(buff+i));

        }
      }

/* FIXME - check if noflags should be made false */

    }
  else
    {
/* main option/setup parse */
/* acquire next option context? */
    code = gulp_is_option(*buff);
    if (code != -1)
      context = code;

/* pass everything we don't recognize to the extra string */
if (code == -1)
  {
/* FIXME - should have finished all contextual stuff */
  if (context != -1)
    {
#if DEBUG_READ_GULP
printf("Warning: context = %d\n", context);
printf(" > %s", scan_cur_line(scan));
#endif
    }
  g_string_sprintfa(ex_buff, "%s", scan_cur_line(scan));
  }

/*
printf("code: %d, context: %d\n", code, context);
*/
/* parse current option */
      switch (code)
        {
        case GULP_COSMO_SHAPE:
          if (num_tokens > 1)
            {
            if (g_strncasecmp(*(buff+1), "dode", 4) == 0)
              gulp.cosmo_shape = 1;
            else
              gulp.cosmo_shape = 0;
            }
          break;
        case GULP_COSMO_POINTS:
          if (num_tokens > 1)
            gulp.cosmo_points = str_to_float(*(buff+1));
          break;
        case GULP_COSMO_SEGMENTS:
          if (num_tokens > 1)
            gulp.cosmo_segments = str_to_float(*(buff+1));
          break;
        case GULP_COSMO_SMOOTHING:
          if (num_tokens > 1)
            gulp.cosmo_smoothing = str_to_float(*(buff+1));
          break;
        case GULP_COSMO_SOLVENT_EPSILON:
          if (num_tokens > 1)
            gulp.cosmo_solvent_epsilon = str_to_float(*(buff+1));
          break;
        case GULP_COSMO_SOLVENT_RADIUS:
          if (num_tokens > 1)
            gulp.cosmo_solvent_radius = str_to_float(*(buff+1));
          if (num_tokens > 2)
            gulp.cosmo_solvent_delta = str_to_float(*(buff+2));
          break;
        case GULP_COSMO_SOLVENT_RMAX:
          if (num_tokens > 1)
            gulp.cosmo_solvent_rmax = str_to_float(*(buff+1));
          break;

        case ENDFORCE:
          g_string_sprintfa(ex_buff, "%s", scan_cur_line(scan));
          break;

        case PRINT:
/* skip, don't store the print option as gdis always appends it */
          if (num_tokens < 2)
            {
            while (!scan_complete(scan))
              {
              g_strfreev(buff);
              buff = scan_get_tokens(scan, &num_tokens);
              if (num_tokens)
                break;
              }
            }
          break;

        case MAXCYC:
          if (num_tokens > 1)
            maxcyc = str_to_float(*(buff+1));
          else
            {
            while (!scan_complete(scan))
              {
              g_strfreev(buff);
              buff = scan_get_tokens(scan, &num_tokens);
              if (num_tokens)
                {
                maxcyc = str_to_float(*buff);
                break;
                }
              }
            }
          break;

        case SWITCH_ON:
          for (j=1 ; j<num_tokens ; j++)
            {
            code = gulp_is_option(*(buff+j));
            switch (code)
              {
              case BFGS_OPT:
              case CONJ_OPT:
              case RFO_OPT:
                optimiser2 = code;
                break;

              case GNORM:
              case CYCLE:
                switch_type = code;
                break;
 
              default:
                if (str_is_float(*(buff+j)))
                  switch_value = str_to_float(*(buff+j));
              }
            }
          break;

/* NB: these keywords can initiate a new model (ie if repeat occurs) */
        case PFRAC:
        case SFRAC:
        case FRAC:
        case CART:
          if (num_tokens > 2)
            {
            if (g_ascii_strncasecmp("region",*(buff+1),6) == 0)
              region = str_to_float(*(buff+2)) - 1;
            model->region_empty[region] = FALSE; 
            if (region == REGION1A)
              model = gulp_change_model(GULP_NEW_COORDS, model);
            }
          else
            model = gulp_change_model(GULP_NEW_COORDS, model);
          nc = ns = 0;
          break;

/* another keyword that can initiate a new model (on repeat) */
        case NAME:
          model = gulp_change_model(GULP_NEW_NAME, model);
          break;

/* other keywords that can initiate a new model (on repeat) */
        case CELL:
        case POLYMER_CELL:
        case POLYMER_VECTOR:
        case SURFACE_CELL:
        case SURFACE_VECTORS:
        case LATTICE_VECTORS:
          model = gulp_change_model(GULP_NEW_LATTICE, model);
          break;
        }
    }

/* main option/data parse */
/* NB: after each case - set context = -1 if parsing for that option is completed */
    switch (context)
      {
/* pass through everything in these blocks */
      case TITLE:
      case WEIGHT:
      case GULP_OBSERVABLES:
      case GULP_VARIABLES:
        g_string_sprintfa(ex_buff, "%s", scan_cur_line(scan));
        while (!scan_complete(scan))
          {
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          if (num_tokens)
            {
/* check first item only */
            code = gulp_is_option(*(buff));
            if (code == -1)
              g_string_sprintfa(ex_buff, "%s", scan_cur_line(scan));
            else
              {
/* terminate when we hit a new option */
/* NB: if option was not an "end" - buffer it for the next parse cycle */
              if (code != END)
                {
                scan_put_line(scan);
                }
              else
                g_string_sprintfa(ex_buff, "%s", scan_cur_line(scan));
              break;
              }
            }
          }
        break;

/* hack for getting an associated .trg file */
      case OUTPUT:
        switch (num_tokens)
          {
          case 4:
            if (g_ascii_strncasecmp(*(buff+1),"mov",3) == 0)
              {
              if (model->gulp.mov_file);
                g_free(model->gulp.mov_file);
              model->gulp.mov_file = g_strdup(*(buff+3));
              }
            break;
          case 3:
            if (g_ascii_strncasecmp(*(buff+1),"traj",4) == 0)
              {
              g_free(model->gulp.trj_file);
              model->gulp.trj_file = g_strdup(*(buff+2));
              model->animation = TRUE;
              }
            break;
          }
        context = -1;
        break;

        case ENSEMBLE:
          if (num_tokens > 1)
            {
            if (g_ascii_strncasecmp(*(buff+1),"nve",3) == 0)
              model->gulp.ensemble = NVE;
            if (g_ascii_strncasecmp(*(buff+1),"nvt",3) == 0)
              model->gulp.ensemble = NVT;
            if (g_ascii_strncasecmp(*(buff+1),"npt",3) == 0)
              model->gulp.ensemble = NPT;
            }
          context = -1;
          break;

        case TEMPERATURE:
          if (num_tokens > 1)
            {
            g_free(model->gulp.temperature);
            model->gulp.temperature = g_strjoinv(" ", buff+1);
            }
          else
            {
/* get next (non empty) line */
            g_strfreev(buff);
            buff = scan_get_tokens(scan, &num_tokens);
            if (buff)
              {
              g_free(model->gulp.temperature);
              model->gulp.temperature = g_strjoinv(" ", buff);
              }
            }
          context = -1;
          break;

        case PRESSURE:
          if (num_tokens > 1)
            {
            g_free(model->gulp.pressure);
            model->gulp.pressure = g_strjoinv(" ", buff+1);
            }
          else
            {
/* get next (non empty) line */
            g_strfreev(buff);
            buff = scan_get_tokens(scan, &num_tokens);
            if (buff)
              {
              g_free(model->gulp.pressure);
              model->gulp.pressure = g_strjoinv(" ", buff);
              }
            }
          context = -1;
          break;

       case TIMESTEP:
         g_free(model->gulp.timestep);
         model->gulp.timestep = g_strdup(get_token_pos(scan_cur_line(scan), 1));
         g_strstrip(model->gulp.timestep);
         break;
       case EQUILIBRATION:
         g_free(model->gulp.equilibration);
         model->gulp.equilibration = g_strdup(get_token_pos(scan_cur_line(scan), 1));
         g_strstrip(model->gulp.equilibration);
         break;
       case PRODUCTION:
         g_free(model->gulp.production);
         model->gulp.production = g_strdup(get_token_pos(scan_cur_line(scan), 1));
         g_strstrip(model->gulp.production);
         break;
       case SAMPLE:
         g_free(model->gulp.sample);
         model->gulp.sample = g_strdup(get_token_pos(scan_cur_line(scan), 1));
         g_strstrip(model->gulp.sample);
         break;
       case WRITE:
         g_free(model->gulp.write);
         model->gulp.write = g_strdup(get_token_pos(scan_cur_line(scan), 1));
         g_strstrip(model->gulp.write);
         break;

        case SUPER_CELL:
          if (num_tokens > 3)
            {
            model->gulp.super[0] = str_to_float(*(buff+1));
            model->gulp.super[1] = str_to_float(*(buff+2));
            model->gulp.super[2] = str_to_float(*(buff+3));
            }
          break;

        case NAME:
          if (num_tokens > 1)
            {
            g_free(model->basename);
            model->basename = g_strjoinv(" ", buff+1);
            g_strstrip(model->basename);
            strcpy(model->filename, model->basename);
            }
          break;

        case GULP_LIBRARY:
          if (num_tokens > 1)
            {
            g_free(model->gulp.libfile);
            model->gulp.libfile = g_strdup(*(buff+1));
            }
          break;

        case SPACE:
          if (num_tokens > 1)
            {
            line = scan_cur_line(scan);
            tmp = get_token_pos(line, 1);
            }
          else
            tmp = scan_get_line(scan);
 
/* process for spacgroup symbol or number */
          gen_flag=0;
          for (i=0 ; i<strlen(tmp) ; i++)
            {
/* any alphabetic chars => space *name* */
            if (g_ascii_isalpha(*(tmp+i)))
              {
              model->sginfo.spacename = g_strdup(tmp);
/* indicate that name should used in lookup */
              model->sginfo.spacenum = -1;
              gen_flag++;
              break;
              }
            }
/* no alphabetic chars, assume space group number */
          if (!gen_flag)
            model->sginfo.spacenum = str_to_float(tmp);
          context = -1;
          break;

        case ORIGIN:
          switch (num_tokens)
            {
            case 2:
              model->sginfo.cellchoice = str_to_float(*(buff+1));
              break;
            default:
              printf("Warning: unsupported origin specification.\n");
            }
          break;

        case POLYMER_CELL:
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          if (num_tokens)
            {
            model->pbc[0] = str_to_float(*buff);
            model->periodic = 1;
            }
          break;

        case POLYMER_VECTOR:
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          if (num_tokens)
            {
            model->latmat[0] = str_to_float(*(buff+0)); 
            model->latmat[3] = 0.0;
            model->latmat[6] = 0.0;
            model->construct_pbc = TRUE;
            model->periodic = 1;
            matrix_lattice_init(model);
            }
          break;

        case SURFACE_CELL:
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          if (num_tokens > 2)
            {
            model->pbc[0] = str_to_float(*(buff+0));
            model->pbc[1] = str_to_float(*(buff+1));
            model->pbc[5] = str_to_float(*(buff+2));
            model->pbc[5] *= D2R;
            model->periodic = 2;
            }
          break;

        case SURFACE_VECTORS:
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          if (num_tokens > 1)
            {
/* FIXME - a is not necessarily colinear with x */
            vec1[0] = str_to_float(*(buff+0));
            vec2[0] = str_to_float(*(buff+1));
            g_strfreev(buff);
            buff = scan_get_tokens(scan, &num_tokens);
            if (num_tokens > 1)
              {
              vec1[1] = str_to_float(*(buff+0));
              vec2[1] = str_to_float(*(buff+1));
/* NB: gdis wants transposed gulp matrix */
              VEC3SET(&(model->latmat[0]), vec1[0], vec1[1], 0.0);
              VEC3SET(&(model->latmat[3]), vec2[0], vec2[1], 0.0);
              VEC3SET(&(model->latmat[6]), 0.0, 0.0, 1.0);
              model->construct_pbc = TRUE;
              model->periodic = 2;
              matrix_lattice_init(model);
              }
            }
          break;

        case CELL:
/* get next line if insufficient tokens */
          if (num_tokens < 2)
            {
            g_strfreev(buff);
            buff = scan_get_tokens(scan, &num_tokens);
            i=0;
            }
          else
            i=1;
/* parse for cell data */
          if (num_tokens > 5)
            {
            model->pbc[0] = str_to_float(*(buff+i+0));
            model->pbc[1] = str_to_float(*(buff+i+1));
            model->pbc[2] = str_to_float(*(buff+i+2));
            model->pbc[3] = str_to_float(*(buff+i+3));
            model->pbc[4] = str_to_float(*(buff+i+4));
            model->pbc[5] = str_to_float(*(buff+i+5));
/* hack to determine if 2D or 3D */
/* FIXME - possible for pbc[3] or pbc[4] to be 0 instead */
            if (fabs(model->pbc[2]) < FRACTION_TOLERANCE)
              model->periodic = 2;
            else
              model->periodic = 3;
            model->pbc[3] *= D2R;
            model->pbc[4] *= D2R;
            model->pbc[5] *= D2R;
            }
          break;

        case LATTICE_VECTORS:
          VEC3SET(vec1, 1.0, 0.0, 0.0);
          VEC3SET(vec2, 0.0, 1.0, 0.0);
          VEC3SET(vec3, 0.0, 0.0, 1.0);
          model->periodic = 0;
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          if (num_tokens > 2)
            {
            vec1[0] = str_to_float(*(buff+0));
            vec2[0] = str_to_float(*(buff+1));
            vec3[0] = str_to_float(*(buff+2));
            model->periodic++;
            }
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          if (num_tokens > 2)
            {
            vec1[1] = str_to_float(*(buff+0));
            vec2[1] = str_to_float(*(buff+1));
            vec3[1] = str_to_float(*(buff+2));
            model->periodic++;
            }
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          if (num_tokens > 2)
            {
            vec1[2] = str_to_float(*(buff+0));
            vec2[2] = str_to_float(*(buff+1));
            vec3[2] = str_to_float(*(buff+2));
            model->periodic++;
            }
          if (model->periodic == 3)
            {
/* use the supplied lattice matrix */
            ARR3SET(&(model->latmat[0]), vec1);
            ARR3SET(&(model->latmat[3]), vec2);
            ARR3SET(&(model->latmat[6]), vec3);
            model->construct_pbc = TRUE;
            matrix_lattice_init(model);
            }
          break;

/* NB: errors -> break (not return) so we cope with empty regions */
        case PFRAC:
        case SFRAC:
        case FRAC:
          model->fractional = TRUE;
        case CART:
/* au check */
          gen_flag = 0;
          if (num_tokens > 1)
            if (g_ascii_strncasecmp("au",*(buff+1),2) == 0)
              gen_flag++;

/* get new line */
/* insufficient items => stupid gulp frac/cart \n number \n coords format */
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          if (num_tokens < 4)
            {
            g_strfreev(buff);
            buff = scan_get_tokens(scan, &num_tokens);
            }

/* core, shell indices */
          while (!scan_complete(scan))
            {
            if (gulp_is_option(*buff) > -1)
              {
              scan_put_line(scan);
              break;
              }

/* NEW - translation and growth slice check */
            growth = translate = FALSE;
            switch (*(*(buff+num_tokens-1)))
              {
              case '%':
                growth = TRUE;
                break;

              case 'T':
                translate = TRUE;
                break;
              }

            if (num_tokens < 4)
              break;

/* if an atom type (ie core/shell) is specified - adjust data positions */
            type = 0;
            breathe = FALSE;
            if (g_ascii_strncasecmp(*(buff+1),"core",4) == 0)
              type=1;
            if (g_ascii_strncasecmp(*(buff+1),"shel",4) == 0)
              type=2;
            if (g_ascii_strncasecmp(*(buff+1),"bcor",4) == 0)
              {
              type=1;
              breathe=TRUE;
              }
            if (g_ascii_strncasecmp(*(buff+1),"bshe",4) == 0)
              {
              type=2;
              breathe=TRUE;
              }

/* read appropriately */
            switch (type)
              {
              case 0:
              case 1:
                if (num_tokens > (3+type))
                  {
                  core = new_core(*buff, model);
                  model->cores = g_slist_prepend(model->cores, core);

                  core->x[0] = str_to_float(*(buff+1+type));
                  core->x[1] = str_to_float(*(buff+2+type));
                  core->x[2] = str_to_float(*(buff+3+type));
                  }
                else
                  break;

                core->breathe = breathe;
                core->region = region;
                core->growth = growth;
                core->translate = translate;

                if (breathe)
                  {
                  if (num_tokens > 5)
                    core->radius = str_to_float(*(buff+5));
/* FIXME - what about charge? */
                  }
                else
                  {
                  i = 4+type; 
                  remaining = num_tokens - i;

                  remaining = CLAMP(remaining, 0, 5);

                  if (noflags)
                    remaining = CLAMP(remaining, 0, 2); 

                  switch (remaining)
                    {
                    case 5:
                      core->lookup_charge = FALSE;
                      core->charge = str_to_float(*(buff+i));
                      core->has_sof = TRUE;
                      core->sof = str_to_float(*(buff+i+1));
                      core->flags = g_strdup_printf("%s %s %s", *(buff+i+2), *(buff+i+3), *(buff+i+4));
                      break;

                    case 4:
                      core->lookup_charge = FALSE;
                      core->charge = str_to_float(*(buff+i));
                      i++;
                    case 3:
                      core->flags = g_strdup_printf("%s %s %s", *(buff+i), *(buff+i+1), *(buff+i+2));
                      break;

                    case 2:
                      core->has_sof = TRUE;
                      core->sof = str_to_float(*(buff+i+1));
                    case 1:
                      core->lookup_charge = FALSE;
                      core->charge = str_to_float(*(buff+i));
                      break;
                    }

                  }
                if (gen_flag)
                  {
                  VEC3MUL(core->x, AU2ANG);
                  }
                nc++;
                break;

              case 2:
                if (num_tokens > 4)
                  {
                  shell = new_shell(*buff, model);
                  model->shels = g_slist_prepend(model->shels, shell);

                  shell->x[0] = str_to_float(*(buff+2));
                  shell->x[1] = str_to_float(*(buff+3));
                  shell->x[2] = str_to_float(*(buff+4));
                  }
                else
                  break;

                shell->breathe = breathe;
                shell->region = region;
                shell->translate = translate;

                if (breathe)
                  {
                  if (num_tokens > 5)
                    shell->radius = str_to_float(*(buff+5));
/* FIXME - what about charge? */
                  }
                else
                  {
                  i = 5;
                  remaining = num_tokens - i;
                  remaining = CLAMP(remaining, 0, 5);

                  if (noflags)
                    remaining = CLAMP(remaining, 0, 2); 

                  switch (remaining)
                    {
                    case 5:
                      shell->lookup_charge = FALSE;
                      shell->charge = str_to_float(*(buff+i));
                      shell->has_sof = TRUE;
                      shell->sof = str_to_float(*(buff+i+1));
                      shell->flags = g_strdup_printf("%s %s %s", *(buff+i+2), *(buff+i+3), *(buff+i+4));
                      break;

                    case 4:
                      shell->lookup_charge = FALSE;
                      shell->charge = str_to_float(*(buff+i));
                      i++;
                    case 3:
                      shell->flags = g_strdup_printf("%s %s %s", *(buff+i), *(buff+i+1), *(buff+i+2));
                      break;

                    case 2:
                      shell->has_sof = TRUE;
                      shell->sof = str_to_float(*(buff+i+1));
                    case 1:
                      shell->lookup_charge = FALSE;
                      shell->charge = str_to_float(*(buff+i));
                      break;
                    }
                  }
                if (gen_flag)
                  {
                  VEC3MUL(shell->x, AU2ANG);
                  }
                ns++;
                break;
              }
            g_strfreev(buff);
            buff = scan_get_tokens(scan, &num_tokens);
            }
          break;

        case D_HKL:
          if (num_tokens > 1)
            model->surface.dspacing = str_to_float(*(buff+1));
          context = -1;
          break;

        case SBULK_ENERGY:
          if (num_tokens > 1)
            model->gulp.sbulkenergy = str_to_float(*(buff+1));
          context = -1;
          break;

        case TOTAL_ENERGY:
          if (num_tokens > 1)
            model->gulp.energy = str_to_float(*(buff+1));
          context = -1;
          break;

        case SPECIES:
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          while (!scan_complete(scan))
            {
/* end on an option */
            if (num_tokens)
              {
              if (gulp_is_option(*buff) > -1)
                {
                scan_put_line(scan);
                break;
                }

/* prevent core dump on comment */
              if (**buff == '#')
                break;

/* add new species data if the 1st item an element, otherwise end */
              if (elem_test(*buff))
                g_string_sprintfa(species_buff, "%s", scan_cur_line(scan));
              else
                break;
              }
            g_strfreev(buff);
            buff = scan_get_tokens(scan, &num_tokens);
            }
          context = -1;
          break;

        case ELEMENT:
          g_strfreev(buff);
          buff = scan_get_tokens(scan, &num_tokens);
          while (!scan_complete(scan))
            {
            if (num_tokens)
              {
/* end on an option */
              temp_code = gulp_is_option(*buff);
              if (temp_code != COVALENT && temp_code != IONIC && temp_code != VDW)
                {
                scan_put_line(scan);
                break;
                }
              else
                g_string_sprintfa(elem_buff, "%s", scan_cur_line(scan));
              }
            g_strfreev(buff);
            buff = scan_get_tokens(scan, &num_tokens);
            }
          context = -1;
          break;

        case KPOINTS:
/* get all subsequent non-empty & numeric lines */
          while (!scan_complete(scan))
            {
            g_strfreev(buff);
            buff = scan_get_tokens(scan, &num_tokens);
            if (num_tokens)
              {
              if (str_is_float(*buff))
                g_string_sprintfa(kpoints_buff, "%s", scan_cur_line(scan));
              }
            else
              break;
            }
          context = -1;
          break;

        case DUMP_FILE:
          if (num_tokens > 1)
            {
/* NEW - grab everything except the dump option itself */
            g_free(model->gulp.dump_file);
            line = scan_cur_line(scan);
            model->gulp.dump_file = g_strdup(get_token_pos(line,1));
            strip_extra(model->gulp.dump_file);
            }
          context = -1;
          break;

        case IGNORE:
/* read until hit an ERONGI */
          while (!scan_complete(scan))
            {
            line = scan_get_line(scan);
            if (gulp_is_option(line) == ERONGI)
              break;
            }
          context = -1;
          break;

/* special case - only one line expected */
        case GULP_SINGLE_LINE_POTENTIAL:
          g_string_sprintfa(ff_buff, "%s", scan_cur_line(scan));
          context = -1;
          break;

/* multi-line potentials */
	case GULP_POTENTIAL:
/* store all potential lines */
          g_string_sprintfa(ff_buff, "%s", scan_cur_line(scan));
          while (!scan_complete(scan))
            {
            line = scan_get_line(scan);

/* ignore (but keep reading) empty/blank lines */
/* NB - keep both tests below (NULL and \n tests) */
            if (!line)
              continue;
            if (strlen(line) < 2)
              continue;
/* terminate on option */
            if (gulp_is_option(line) > -1)
              {
              context = -1;
              scan_put_line(scan);
              break;
              }
/* NEW - terminate if some chars of 1st token are alphabetic and not an element */
            g_strfreev(buff);
            buff = tokenize(line, &num_tokens);
            if (num_tokens)
              {
              tmp = NULL;
              for (i=strlen(*buff) ; i-- ; )
                {
                if (g_ascii_isalpha((*buff)[i]))
                  {
                  tmp = *buff;
                  break;
                  }
                }
              if (tmp)
                {
/* NEW - allow GULP wildcard element */
                if (g_ascii_strncasecmp(tmp, "x\0", 2) != 0)
                  {
                  if (!elem_symbol_test(tmp))
                    {
                    context = -1;
                    scan_put_line(scan);
                    break;
                    }
                  }
                } 
              }
/* store the string */
            g_string_sprintfa(ff_buff, "%s", line);
            }
          break;
        }

/* get next non-empty line */
/* TODO - encapsulate in function call? (buff & line) */
  g_strfreev(buff);
  buff = scan_get_tokens(scan, &num_tokens);
  }
g_strfreev(buff);

/* TODO - if keywords found > 0, but doesn't match num tokens */
/* then we've found an unrecognized keyword -> pass through?/print warning? */
/*
printf("Keywords processed: %d/%d\n", keywords_found, keywords_expected);
*/

/* NEW */
if (gulp.solvation_model != GULP_SOLVATION_NONE)
  if (!gulp_shape_indices_compute(&gulp))
    printf("WARNING: input points may be invalid.\n");

/* setup newly loaded models */
for (item=sysenv.mal ; item ; item=g_slist_next(item))
  {
  model = item->data;

  if (model->id == -1)
    {
    model->id = GULP;
    strcpy(model->filename, filename);

/* NEW - unrecognized keyword transfer */
    if (keywords_found != keywords_expected)
      {
      g_free(model->gulp.extra_keywords);
      model->gulp.extra_keywords = g_strdup(key_buff->str);
      }

/* no depth info for loaded surfaces, so pretend they're MARVIN style */
    if (model->periodic == 2)
      model->surface.true_cell = FALSE;

/* NEW */
model->gulp.solvation_model = gulp.solvation_model;
model->gulp.cosmo_shape = gulp.cosmo_shape;
model->gulp.cosmo_shape_index[0] = gulp.cosmo_shape_index[0];
model->gulp.cosmo_shape_index[1] = gulp.cosmo_shape_index[1];
model->gulp.cosmo_points = gulp.cosmo_points;
model->gulp.cosmo_segments = gulp.cosmo_segments;
model->gulp.cosmo_smoothing = gulp.cosmo_smoothing;
model->gulp.cosmo_solvent_epsilon = gulp.cosmo_solvent_epsilon;
model->gulp.cosmo_solvent_radius = gulp.cosmo_solvent_radius;
model->gulp.cosmo_solvent_delta = gulp.cosmo_solvent_delta;
model->gulp.cosmo_solvent_rmax = gulp.cosmo_solvent_rmax;


/* transfer universal gulp data */
/* FIXME - deprec - use gulp_pak instead (see above cosmo stuff) */
    model->gulp.run = run;
    model->gulp.method = method;
    model->gulp.optimiser = optimiser;
    if (optimiser2 > 0)
      {
      model->gulp.optimiser2 = optimiser2;
      model->gulp.switch_type = switch_type;
      model->gulp.switch_value = switch_value;
      }
    model->gulp.unit_hessian = unit;
    model->gulp.coulomb = coulomb;
    model->gulp.maxcyc = maxcyc;
    model->gulp.qeq = qeq;
    model->gulp.free = free;
    model->gulp.zsisa = zsisa;
    model->gulp.compare = compare;
    model->gulp.nosym = nosym;
    model->gulp.fix = fix;
    model->gulp.phonon = phonon;
    model->gulp.eigen = eigen;

    model->gulp.kpoints = g_strdup(kpoints_buff->str);
    model->gulp.potentials = g_strdup(ff_buff->str);

/* duplicate species data */
    model->gulp.species = g_strdup(species_buff->str);

    buff = tokenize(model->gulp.species, &num_tokens);
    i = 0;
    while (i < num_tokens)
      {
      type = elem_test(*(buff+i));
      j = i;

      get_elem_data(type, &elem_data, model);
      if (type && ++i < num_tokens)
        {
        if (g_ascii_strncasecmp(*(buff+i), "shel", 4) == 0)
          {
          if (++i < num_tokens)
            {
            elem_data.shell_charge = str_to_float(*(buff+i));
            for (slist=model->shels ; slist ; slist=g_slist_next(slist))
              {
              shell = slist->data;
              if (shel_match(*(buff+j), shell))
                {
/* NEW - only assign species data if explict charges were absent */
                if (shell->lookup_charge)
                  {
                  shell->charge = elem_data.shell_charge;
                  shell->lookup_charge = FALSE;
                  }
                }
              }
            }
          }
        else
          {
          if (g_ascii_strncasecmp(*(buff+i), "core", 4) == 0)
            {
            if (++i < num_tokens)
              elem_data.charge = str_to_float(*(buff+i));
            else
              elem_data.charge = 0.0;
            }
          else
            elem_data.charge = str_to_float(*(buff+i));
          for (clist=model->cores ; clist ; clist=g_slist_next(clist))
            {
            core = clist->data;
            if (core_match(*(buff+j), core))
              {
/* NEW - only assign species data if explict charges were absent */
              if (core->lookup_charge)
                {
                core->charge = elem_data.charge;
                core->lookup_charge = FALSE;
                }
              }
            }
          }
        } 
      i++;
      }
    g_strfreev(buff);

/* process element data */
    model->gulp.elements = g_strdup(elem_buff->str);
    buff = tokenize(model->gulp.elements, &num_tokens);
    i = 0;
    type = -1;
    while (i < num_tokens)
      {
      switch (type)
        {
        case COVALENT:
          code = elem_test(*(buff+i));
          i++;
          if (!get_elem_data(code, &elem_data, model) && i<num_tokens)
            {
            elem_data.cova = str_to_float(*(buff+i));
            put_elem_data(&elem_data, model);
            }
          type = -1;
          break;

        case VDW:
          code = elem_test(*(buff+i));
          i++;
          if (!get_elem_data(code, &elem_data, model) && i<num_tokens)
            {
            elem_data.vdw = str_to_float(*(buff+i));
            put_elem_data(&elem_data, model);
            }
          type = -1;
          break;

        default:
          if (g_ascii_strncasecmp(*(buff+i), "cova", 4) == 0)
            type = COVALENT;
          if (g_ascii_strncasecmp(*(buff+i), "vdw", 3) == 0)
            type = VDW;
        }
      i++;
      }
    g_strfreev(buff);

/* extra (unprocessed lines) */
    if (strlen(ex_buff->str))
      model->gulp.extra = g_strdup(ex_buff->str);

/* NEW - store initial frame coord type, as trajectory animation will overwrite */
    model->gulp.orig_fractional = model->fractional;
    model->gulp.orig_construct_pbc = model->construct_pbc;

/* ensure the ordering is correct, as GULP trajectory files depend on it */
    model->cores = g_slist_reverse(model->cores);
    model->shels = g_slist_reverse(model->shels);

/* display init */
    model_prep(model);
    }
  }


/* cleanup */
g_string_free(key_buff, TRUE);
g_string_free(elem_buff, TRUE);
g_string_free(ex_buff, TRUE);
g_string_free(ff_buff, TRUE);
g_string_free(species_buff, TRUE);
g_string_free(kpoints_buff, TRUE);

scan_free(scan);
return(0);
}

/*******************************************************************/
/* primitive for adding two strings to generate the property value */
/*******************************************************************/
void property_add_with_units(gint rank,
                             gchar *key, gchar *value, gchar *units,
                             struct model_pak *model)
{
gchar *text;

text = g_strjoin(" ", value, units, NULL);
property_add_ranked(rank, key, text, model);
g_free(text);
}

/*****************************/
/* generic gulp output parse */
/*****************************/
/* NB: new model is not made, as only energy etc. values */
/* are of interest if we're parsing the results of a gulp run */
#define DEBUG_READ_GULP_OUTPUT 0
gint read_gulp_output(gchar *filename, struct model_pak *data)
{
gint i, j, m, n, flag, model, region, primitive=FALSE, comp=0;
gint status, count=0, dflag=0, fflag=0, pflag=0, read_pass=0, num_tokens;
gchar line[LINELEN], **buff, *text, coord;
gdouble *value;
GSList *list, *flist=NULL, *ir_list=NULL, *raman_list=NULL, *model_list=NULL;
GHashTable *types=NULL;
struct model_pak *temp;
struct core_pak *core;
struct shel_pak *shell;
FILE *fp;

/* checks */
g_assert(data != NULL);
g_assert(filename != NULL);

fp = fopen(filename, "rt");
if (!fp)
  return(1);

/* init */
model = data->number;
data->region_empty[REGION1A] = FALSE;
data->gulp.esurf[0] = 0.0;
data->gulp.esurf[1] = 0.0;
data->gulp.eatt[0] = 0.0;
data->gulp.eatt[1] = 0.0;

types = g_hash_table_new_full(g_str_hash, g_str_equal, g_free, g_free);

#if DEBUG_READ_GULP_OUTPUT
printf("Grafted: %d\n", data->grafted);
#endif

/* init for new model (ie not on the tree) */
if (!data->grafted)
  {
  data->id = GULPOUT;
  g_free(data->basename);
  data->basename = strdup_basename(filename);
  strcpy(data->filename, filename);
  model_list = g_slist_prepend(model_list, data);
  }

/* look for interesting output */
flag=status=0;
while(!fgetline(fp,line))
  {
  buff = tokenize(line, &num_tokens);

/* GULP trapping */
  if (g_ascii_strncasecmp("!! ERROR",line,8) == 0)
    {
    data->error_file_read = g_string_append(data->error_file_read, line);
    g_strfreev(buff);
    status = 2;
    break;
    }
  if (g_ascii_strncasecmp("  **** Unit cell is not charge neutral",line,38) == 0)
    {
    data->error_file_read = g_string_append(data->error_file_read, line);
    if (!fgetline(fp,line))
      data->error_file_read = g_string_append(data->error_file_read, line);
    g_strfreev(buff);
    status = 2;
    break;
    }
/*
  if (g_ascii_strncasecmp("  **** Warning", line, 14) == 0)
    gui_text_show(WARNING, line);
*/

/* periodicity */
  if (g_ascii_strncasecmp("  Dimensionality", line, 16) == 0)
    {
    if (num_tokens > 2)
      data->periodic = str_to_float(*(buff+2));
    }

/* space group */
  if (g_ascii_strncasecmp("  Space group ", line, 14) == 0)
    {
    if (!data->grafted)
      {
      for (i=strlen(line) ; i-- ; )
        {
        if (line[i] == ':')
          {
          data->sginfo.spacenum = -1;
          data->sginfo.spacename = g_strstrip(g_strdup(&line[i+1]));

/*
          if (data->sginfo.spacename[0] == 'R')
            primitive = TRUE;
      if (g_strrstr(data->sginfo.spacename, ":R") )
        {
printf("1 prim\n");
        primitive = TRUE;
        }
*/

#if DEBUG_READ_GULP_OUTPUT
printf("space group: %s\n", data->sginfo.spacename);
#endif
          }
        }
      }
    else
      {
/* NEW - check if we should be looking for the primitve cell */
/*
      if (g_strrstr(data->sginfo.spacename, ":R") )
        {
printf("2 prim\n");
        primitive = TRUE;
        }
*/

      }
    }

/* NEW - process GULP's library symbols as FF types */
  if (g_ascii_strncasecmp("  Species output", line, 16) == 0)
    {
    for (i=4 ; i-- ; )
      {
      g_strfreev(buff);
      buff = get_tokenized_line(fp, &num_tokens);
      }

    while (*buff)
      {
      g_strfreev(buff);
      buff = get_tokenized_line(fp, &num_tokens);

      if (elem_symbol_test(*buff))
        {
        if (num_tokens > 8)
          {
          g_hash_table_insert(types, g_strdup(*buff), g_strdup(*(buff+8)));
          }
        }
      else
        break; 
      }
    }

/* create new model if the calculation was a defect calc */
  if (g_ascii_strncasecmp("*  Defect calculation for configuration", line, 39) == 0)
    {
/* create new model */
    temp = model_new();
    temp->id = GULPOUT;
    g_free(temp->basename);
    temp->basename = g_strdup_printf("%s_defect_r1", data->basename);
/* this is now our new model */
    data = temp;
    model_list = g_slist_prepend(model_list, data);
    }

/* create new model on multiple instances of output data */
/* introducted to cope with output from the translate option */
  if (g_ascii_strncasecmp("  Components of energy", line, 22) == 0)
    {
    if (comp && !(comp % 2))
      {
/* create new model */
      temp = model_new();
      temp->id = GULPOUT;
/* pass model data through */
      temp->fractional = data->fractional;
      temp->periodic = data->periodic;
      temp->sginfo.spacename = g_strdup(data->sginfo.spacename);
      temp->sginfo.spacenum = -1;
      temp->construct_pbc = data->construct_pbc;
      if (data->construct_pbc)
        memcpy(temp->latmat, data->latmat, 9*sizeof(gdouble));
      else
        memcpy(temp->pbc, data->pbc, 6*sizeof(gdouble));
/* this is now our new model */
      data = temp;
      model_list = g_slist_prepend(model_list, data);
      }
    comp++;
    }

/* minimized energy search */
  if (g_ascii_strncasecmp("  Final energy =",line,16) == 0)
    {
    if (num_tokens > 3)
      {
/* record total energy */
      data->gulp.energy = str_to_float(*(buff+3));
      flag++;
      }
    }

/* surface energy search */
  if (g_ascii_strncasecmp("  Surface energy ",line,17) == 0)
    {
    if (num_tokens > 6)
      {
      if (data->gulp.esurf[0] == 0.0)
        {
        data->gulp.esurf[0] = str_to_float(*(buff+5));
        property_add_with_units(3, "Esurf (u)", *(buff+5), *(buff+6), data);
        }
      else
        {
        data->gulp.esurf[1] = str_to_float(*(buff+5));
        property_add_with_units(4, "Esurf (r)", *(buff+5), *(buff+6), data);
        }
      g_free(data->gulp.esurf_units);
      data->gulp.esurf_units = g_strdup(*(buff+6));
      }
    }

/* attachment energy search 1 */
  if (g_ascii_strncasecmp("  Attachment energy ",line,20) == 0)
    {
    if (num_tokens > 4)
      {
      if (data->gulp.eatt[0] == 0.0)
        {
        data->gulp.eatt[0] = str_to_float(*(buff+3));
        property_add_with_units(5, "Eatt (u)", *(buff+3), *(buff+4), data);
        }
      else
        {
        data->gulp.eatt[1] = str_to_float(*(buff+3));
        property_add_with_units(6, "Eatt (r)", *(buff+3), *(buff+4), data);
        }
      g_free(data->gulp.eatt_units);
      data->gulp.eatt_units = g_strdup(*(buff+4));
      }
    }

/* attachment energy search 2 */
  if (g_ascii_strncasecmp("  Attachment energy/",line,20) == 0)
    {
    if (num_tokens > 3)
      {
/* NB: assumes this match occurs just after the previous match */
      if (data->gulp.eatt[1] == 0.0)
        {
        data->gulp.eatt[0] = str_to_float(*(buff+3));
        property_add_with_units(5, "Eatt (u)", *(buff+3), "eV/unit", data);
        }
      else
        {
        data->gulp.eatt[1] = str_to_float(*(buff+3));
        property_add_with_units(6, "Eatt (r)", *(buff+3), "eV/unit", data);
        }
      g_free(data->gulp.eatt_units);
      data->gulp.eatt_units = g_strdup("eV/unit");
      }
    }

/* search for gnorm */
  if (g_ascii_strncasecmp("  Final Gnorm  =",line,16) == 0)
    {
    if (num_tokens > 3)
      {
      data->gulp.gnorm = str_to_float(*(buff+3));
      property_add_ranked(9, "gnorm", *(buff+3), data);
      flag++;
      }
    }

/* single point energy search */
  if (g_ascii_strncasecmp("  Total lattice energy",line,22) == 0)
    {
/* NB: assumes we have the correct cell parameters */
    primitive = space_primitive_cell(data);

/* if nothing on the rest of this line, then should be on the next */
    switch (num_tokens)
      {
      case 4:
      while(!fgetline(fp,line))
        {
        if (primitive)
          {
          if (g_ascii_strncasecmp("    Primitive unit cell",line,23) == 0)
            {
            g_strfreev(buff);
            buff = tokenize(line, &num_tokens);
            if (num_tokens > 5)
              data->gulp.energy = str_to_float(*(buff+4));
            break;
            }
          }
        else
          {
          if (g_ascii_strncasecmp("    Non-primitive unit cell",line,27) == 0)
            {
            g_strfreev(buff);
            buff = tokenize(line, &num_tokens);
            if (num_tokens > 5)
              data->gulp.energy = str_to_float(*(buff+4));
            break;
            }
          }
        }
      break;

      case 6:
        if (g_ascii_strncasecmp("eV",*(buff+5),2) == 0)
          data->gulp.energy = str_to_float(*(buff+4));
      break;

      }
    }

/* single point free energy search */
  if (g_ascii_strncasecmp("  Initial surface dipole ",line,25) == 0)
    {
    if (num_tokens > 4)
      data->gulp.sdipole = str_to_float(*(buff+4));
    }

/* single point free energy search */
  if (g_ascii_strncasecmp("  Total free energy",line,19) == 0)
    {
    if (num_tokens > 5)
      if (g_ascii_strncasecmp("eV",*(buff+5),2) == 0)
        data->gulp.energy = str_to_float(*(buff+4));
    }

/* NEW - vibrational info search */
  if (g_ascii_strncasecmp(" Frequencies (cm-1) and Eigenvectors",line,36) == 0)
    fflag++;
  if (g_ascii_strncasecmp("  Zero point energy",line,19) == 0)
    fflag=0;

  if (fflag && g_ascii_strncasecmp(" Frequency",line,10) == 0)
    {
/* space group processing here, since we want to attach */
/* the vibration lists to ALL atoms in the full cell, */
/* not just the asymmetric ones. */
    if (!data->grafted && !pflag)
      {
      model_prep(data);
      pflag++;
      }

    for (i=1 ; i<num_tokens ; i++)
      flist = g_slist_prepend(flist, g_strdup(*(buff+i)));

/* skip to data */
    while (!fgetline(fp,line))
      {
      g_strfreev(buff);
      buff = tokenize(line, &n);

/* NEW - store IR intensities */
      if (g_ascii_strncasecmp(" IR Intensity", line, 13) == 0)
        {
        i = 2;
        for (i=2 ; i<n ; i++)
          ir_list = g_slist_prepend(ir_list, g_strdup(*(buff+i)));
        }

/* NEW - store Raman intensities */
      if (g_ascii_strncasecmp(" Raman", line, 6) == 0)
        {
        i = 2;
        for (i=2 ; i<n ; i++)
          raman_list = g_slist_prepend(raman_list, g_strdup(*(buff+i)));
        }

      if (n > 1)
        coord = **(buff+1);
      else
        coord = '?';

      if (coord == 'x')
        break;
      }

/* read the eigenvectors in */
    m=0;
    while (m<3*g_slist_length(data->cores))
      {
      g_strfreev(buff);
      buff = tokenize(line, &n);
/* blank line termination */
      if (!n)
        break;
/* get reference atom number (NB: GULP starts at 1, list starts at 0) */
      i = str_to_float(*buff);
      i--;

/* FIXME - if we're reading in cores from this gulp file - order will be reverse */
/* due to prepending to core list */
      core = g_slist_nth_data(data->cores, i);
      g_assert(core != NULL);

/*
      if (i<0 || i>data->num_atoms-1)
        {
        printf("Bad line: %s\n", line);
        printf("ref atom: %d, num atoms: %d\n", i, data->num_atoms);
        g_assert_not_reached();
        }
*/
      for (j=2 ; j<n ; j++)
        {
/* get the value */
        value = g_malloc(sizeof(gdouble));
        *value = str_to_float(*(buff+j));
/* get reference coordinate */
        coord = **(buff+1);
        switch(coord)
          {
          case 'x':
            core->vibx_list = g_slist_append(core->vibx_list,
                                         (gpointer *) value);
            break;
          case 'y':
            core->viby_list = g_slist_append(core->viby_list,
                                         (gpointer *) value);
            break;
          case 'z':
            core->vibz_list = g_slist_append(core->vibz_list,
                                         (gpointer *) value);
            break;
          default:
            g_assert_not_reached();
          }
        }

/* next line (if available) */
      if (fgetline(fp,line))
        break;
      m++;
      }
    }

/* minimized coords search */
/* TODO - unminimized coords? */
/* added 'of surface' to avoid reading in growth slice */
  dflag=0;
  if (g_ascii_strncasecmp("  Final asymmetric unit coordinates",line,35) == 0
   || g_ascii_strncasecmp("  Final fractional coordinates",line,30) == 0
   || g_ascii_strncasecmp("  Fractional coordinates of asymmetric",line,38) == 0
   || g_ascii_strncasecmp("  Final fractional/Cartesian coordinates", line, 40) == 0
   || g_ascii_strncasecmp("  Mixed fractional/Cartesian coordinates", line, 40) == 0)
    {
    data->fractional = TRUE;
    dflag++;
    }

/* don't read in coords if already on the tree (eg single pt) */
  if (dflag && !data->grafted)
    {
/* enforce empty core/shell lists */
    free_core_list(data);
/* skip the two (minimum number) ----- dividers */
    region = REGION1A;
    count=0;
    while (count < 2)
      {
      if (fgetline(fp,line))
        return(4);
      if (g_ascii_strncasecmp("-------", line, 7) == 0)
        count++;
      }

/* loop until we run out of atomic data lines */
    i=j=n=0;
    flag=0;
/* read as many atoms as we can */
    while(!fgetline(fp,line))
      {
      g_strfreev(buff);
      buff = tokenize(line, &num_tokens);

/* TODO - dont exit if ----- (could be region divider) */
/*
      if (g_ascii_strncasecmp("-------", line, 7) == 0)
        continue;
*/

/* blank line termination */
if (!num_tokens)
  break;

/* determine region */
      if (g_ascii_strncasecmp("  Region ", line, 9) == 0)
        {
        if (num_tokens > 1)
        switch(*(*(buff+1)))
          {
          case '1':
#if DEBUG_READ_GULP_OUTPUT
printf("Labelling region 1...\n");
#endif
            region = REGION1A;
            break;
          case '2':
#if DEBUG_READ_GULP_OUTPUT
printf("Labelling region 2...\n");
#endif
            region = REGION2A;
            break;
          default:
            printf("WARNING: unknown region specification\n");
            region = REGION1A;
          }
        data->region_empty[region] = FALSE;
        continue;
        }

/* HACK - GULP sometimes puts *'s after coords (region 1 only?) */
      if (num_tokens > 7)
      if (g_ascii_strncasecmp("*", *(buff+4), 1) == 0)
        {
        g_free(*(buff+4)); 
        *(buff+4) = g_strdup(*(buff+5));
        g_free(*(buff+5)); 
        *(buff+5) = g_strdup(*(buff+7));
        }

/* exit only if insufficient data items (NB: changes b/w sections, min=6) */
      if (num_tokens > 6)
        {
/* read the data in - assume the same order */
        if (g_ascii_strncasecmp("c", *(buff+2), 1) == 0)
          {
          core = new_core(*(buff+1), data);
          data->cores = g_slist_prepend(data->cores, core);

          if (!read_pass)
            core->region = region;
          core->x[0] = str_to_float(*(buff+3));
          core->x[1] = str_to_float(*(buff+4));
          core->x[2] = str_to_float(*(buff+5));
          core->charge = str_to_float(*(buff+6));
          core->lookup_charge = FALSE;
/*
          core->sof = str_to_float(*(buff+7));
*/
          i++;

#if DEBUG_READ_GULP_OUTPUT
printf("coords: %lf %lf %lf\n", core->x[0], core->x[1], core->x[2]);
#endif
          }

        if (g_ascii_strncasecmp("s", *(buff+2), 1) == 0)
          {
          shell = new_shell(*(buff+1), data);
          data->shels = g_slist_prepend(data->shels, shell);

          if (!read_pass)
            shell->region = region;
          shell->x[0] = str_to_float(*(buff+3));
          shell->x[1] = str_to_float(*(buff+4));
          shell->x[2] = str_to_float(*(buff+5));
          shell->charge = str_to_float(*(buff+6));
          shell->lookup_charge = FALSE;

#if DEBUG_READ_GULP_OUTPUT
printf("coords: %lf %lf %lf\n", shell->x[0], shell->x[1], shell->x[2]);
#endif

          j++;
          }
        }

      }

data->cores = g_slist_reverse(data->cores);
data->shels = g_slist_reverse(data->shels);

/* done one pass => read in unrelaxed coords if reading optimization */
    read_pass++;

#if DEBUG_READ_GULP_OUTPUT
printf("Retrieved %d atoms & %d shells (periodic).\n",i,j);
#endif
    }

/* cell parameters - search 1 */
  if (g_ascii_strncasecmp("  Final cell parameters",line,22) == 0 && !data->grafted)
    {
/* skip to data */
    fgetline(fp,line);
    fgetline(fp,line);
/* get cell lengths */
    for (i=0 ; i<6 ; i++)
      {
      if (fgetline(fp,line))
        break;
      g_strfreev(buff);
      buff = tokenize(line, &num_tokens);
      if (num_tokens > 1)
        data->pbc[i] = str_to_float(*(buff+1));
      }
/* convert to radians */
    data->pbc[3] *= D2R;
    data->pbc[4] *= D2R;
    data->pbc[5] *= D2R;
    }

/* cell parameters - search 2 */
  if (g_ascii_strncasecmp("  Non-primitive lattice parameters",line,34) == 0 && !data->grafted)
    {
/* skip to data */
    if (fgetline(fp,line))
      break;
    if (fgetline(fp,line))
      break;

/* get cell lengths */
    g_strfreev(buff);
    buff = tokenize(line, &num_tokens);
    if (num_tokens > 8)
      {
      data->pbc[0] = str_to_float(*(buff+2));
      data->pbc[1] = str_to_float(*(buff+5));
      data->pbc[2] = str_to_float(*(buff+8));
      }   
/* get cell angles */
    if (fgetline(fp,line))
      break;
    g_strfreev(buff);
    buff = tokenize(line, &num_tokens);
    if (num_tokens > 5)
      {
      data->pbc[3] = D2R*str_to_float(*(buff+1));
      data->pbc[4] = D2R*str_to_float(*(buff+3));
      data->pbc[5] = D2R*str_to_float(*(buff+5));
      }
    }

/* cell parameters - search 3 */
  if (g_ascii_strncasecmp("  Final surface cell parameters ",line,32) == 0 && !data->grafted)
    {
    data->periodic = 2;
    data->fractional = TRUE;
    data->pbc[2] = 0.0;
    data->pbc[3] = PI/2.0;
    data->pbc[4] = PI/2.0;

/* skip to 1st line of data */
    fgetline(fp,line);
    fgetline(fp,line);

    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 1)
      data->pbc[0] = str_to_float(*(buff+1));

    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 1)
      data->pbc[1] = str_to_float(*(buff+1));

    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 1)
      data->pbc[5] = D2R*str_to_float(*(buff+1));
    }

/* cell parameters - search 4 (eg a CONV surface calc) */
  if (g_ascii_strncasecmp("  Surface cell parameters ",line, 26) == 0 && !data->grafted)
    {
    data->periodic = 2;
    data->fractional = TRUE;
    data->pbc[2] = 0.0;
    data->pbc[3] = PI/2.0;
    data->pbc[4] = PI/2.0;
/* skip to 1st line of data */
    fgetline(fp,line);

    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 5)
      {
      data->pbc[0] = str_to_float(*(buff+2));
      data->pbc[5] = D2R*str_to_float(*(buff+5));
      }
    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 1)
      data->pbc[1] = str_to_float(*(buff+2));
    }

/* cell parameters - search 5 (unrelaxed) */
  if (g_ascii_strncasecmp("  Cell parameters (Angstroms/Degrees):",line,38) == 0 && !data->grafted)
    {
/* skip blank line */
    fgetline(fp,line);
/* get cell lengths */
    for (i=0 ; i<3 ; i++)
      {
      g_strfreev(buff);
      buff = get_tokenized_line(fp, &num_tokens);

      if (num_tokens > 5)
        {
        data->pbc[i] = str_to_float(*(buff+2));
        data->pbc[i+3] = D2R*str_to_float(*(buff+5));
        }
      }
    }

/* cell parameters - search 6 */
  if (g_ascii_strncasecmp("  Surface Cartesian vectors (Angstroms) :",line,41) == 0 && !data->grafted)
    {
    data->periodic = 2;
    data->construct_pbc = TRUE;
/* skip blank line */
    fgetline(fp,line);

/* get vec 1 */
    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 1)
      {
      data->latmat[0] = str_to_float(*buff);
      data->latmat[3] = str_to_float(*(buff+1));
      data->latmat[6] = 0.0;
      }

/* get vec 2 */
    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 1)
      {
      data->latmat[1] = str_to_float(*buff);
      data->latmat[4] = str_to_float(*(buff+1));
      data->latmat[7] = 0.0;
      }

/* set vec 3 */
    data->latmat[2] = 0.0;
    data->latmat[5] = 0.0;
    data->latmat[8] = 1.0;
    }

/* cell parameters - search 8 */
  if (g_ascii_strncasecmp("  Primitive cell parameters :",line,29) == 0 && !data->grafted)
    {
    data->periodic = 3;
/* skip blank line */
    fgetline(fp,line);

/* get line 1 */
    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 11)
      {
      data->pbc[0] = str_to_float(*(buff+8));
      data->pbc[3] = D2R*str_to_float(*(buff+11));
      }
/* get line 2 */
    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 11)
      {
      data->pbc[1] = str_to_float(*(buff+8));
      data->pbc[4] = D2R*str_to_float(*(buff+11));
      }
/* get line 3 */
    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 11)
      {
      data->pbc[2] = str_to_float(*(buff+8));
      data->pbc[5] = D2R*str_to_float(*(buff+11));
      }
    }

/* cell parameters - search 9 */
  if (g_ascii_strncasecmp("  Polymer cell parameter",line,24) == 0 && !data->grafted)
    {
/* init */
    data->periodic = 1;
    data->construct_pbc = FALSE;
/* blank line skip */
    fgetline(fp,line);

    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 2)
      data->pbc[0] = str_to_float(*(buff+2));
    }

/* cell parameters - search 10 */
  if (g_ascii_strncasecmp("  Final polymer cell parameter",line,30) == 0 && !data->grafted)
    {
/* init */
    data->periodic = 1;
    data->construct_pbc = FALSE;
/* blank line skip */
    fgetline(fp,line);
    fgetline(fp,line);

    g_strfreev(buff);
    buff = get_tokenized_line(fp, &num_tokens);
    if (num_tokens > 2)
      data->pbc[0] = str_to_float(*(buff+1));
    }

/* isolated coords search */
  if ((g_ascii_strncasecmp("  Final cartesian coordinates", line, 29) == 0 
   ||  g_ascii_strncasecmp("  Cartesian coordinates of cluster", line, 34) == 0
   ||  g_ascii_strncasecmp("  Region 1 (absolute coordinates)", line, 33) == 0)
   && !data->grafted)
    {
    data->periodic = 0;
    data->fractional = FALSE;

/* enforce empty core/shell lists */
    free_core_list(data);

    for (i=0 ; i<5 ; i++)
      fgetline(fp,line);
/* loop until we run out of atomic data lines */
    i=j=n=0;
    flag=0;
/* don't read any more than num_atoms -> memory allocation problems! */
    while(!fgetline(fp,line) && !flag)
      {
      g_strfreev(buff);
      buff = tokenize(line, &num_tokens);

/* blank line termination */
if (!num_tokens)
  break;

/* exit if atom number is not consecutive */ 
      m = str_to_float(*buff);
      if (m != ++n)
        flag=1;
      else
        {
/* read the data in - assume the same order */
        if (num_tokens > 6)
          {
          if (g_ascii_strncasecmp("c",*(buff+2),1) == 0)
            {
            core = new_core(*(buff+1), data);
            data->cores = g_slist_prepend(data->cores, core);

            core->x[0] = str_to_float(*(buff+3));
            core->x[1] = str_to_float(*(buff+4));
            core->x[2] = str_to_float(*(buff+5));
            core->charge = str_to_float(*(buff+6));
            core->lookup_charge = FALSE;
            }

          if (g_ascii_strncasecmp("s",*(buff+2),1) == 0)
            {
            shell = new_shell(*(buff+1), data);
            data->shels = g_slist_prepend(data->shels, shell);

            shell->x[0] = str_to_float(*(buff+3));
            shell->x[1] = str_to_float(*(buff+4));
            shell->x[2] = str_to_float(*(buff+5));
            shell->charge = str_to_float(*(buff+6));
            shell->lookup_charge = FALSE;
            }
          }
        }
      }

    data->cores = g_slist_reverse(data->cores);
    data->shels = g_slist_reverse(data->shels);

#if DEBUG_READ_GULP_OUTPUT
printf("Retrieved %d atoms & %d shells (cluster).\n",i,j);
#endif
    }

/* NEW */
  if (g_ascii_strncasecmp("  Electrostatic potential at input sites", line, 40) == 0)
    {
/* skip to data */
    for (i=5 ; i-- ; )
      fgetline(fp,line);

/* clean the list */
    if (data->gulp.epot_vals)
      free_slist(data->gulp.epot_vals);
    data->gulp.epot_vals = NULL;

/* acquire lines */
    while (!fgetline(fp, line))
      {
      if (g_ascii_strncasecmp("----------",line,10) == 0)
        break;

      g_strfreev(buff);
      buff = tokenize(line, &num_tokens);
      if (num_tokens > 4)
        {
/* get potential value */
        value = g_malloc(sizeof(gdouble));
        *value = str_to_float(*(buff+4));
        data->gulp.epot_vals = g_slist_prepend(data->gulp.epot_vals, value);
/* record min & max */
        if (*value > data->gulp.epot_max && data->epot_autoscale)
          data->gulp.epot_max = *value;
        if (*value < data->gulp.epot_min && data->epot_autoscale)
          data->gulp.epot_min = *value;
        }
      }

/* NB: prepend -> reverse to get matching direction with epot_vecs */
    data->gulp.epot_vals = g_slist_reverse(data->gulp.epot_vals); 

/* confirm that we got sufficient points */
#if DEBUG_READ_GULP_OUTPUT
    printf("got: %d, expected: %d\n",
           g_slist_length(data->gulp.epot_vecs),
           g_slist_length(data->gulp.epot_vals));
    printf("min: %f, max: %f\n", data->gulp.epot_min, data->gulp.epot_max);
#endif
    }

  g_strfreev(buff);
  }

#if DEBUG_READ_GULP_OUTPUT
P3VEC("cell lengths ",&data->pbc[0]);
P3VEC("cell angles ",&data->pbc[3]);
#endif

text = g_strdup_printf("%.5f eV", data->gulp.energy);
property_add_ranked(2, "Energy", text, data);
g_free(text);

/* surfaces are always const vol */
if (data->periodic == 2)
  data->gulp.method = CONV;

/* init lists */
if (flist)
  data->phonons = g_slist_reverse(flist);
else
  data->phonons = NULL;
data->num_phonons = g_slist_length(data->phonons);

if (ir_list)
  data->ir_list = g_slist_reverse(ir_list); 
if (raman_list)
  data->raman_list = g_slist_reverse(raman_list); 

#if DEBUG_READ_GULP_OUTPUT
printf("vibrational eigenvectors: %d\n", data->num_phonons);
if (data->num_phonons)
  dump_phonons(data);
#endif

/* initialize the models for display */
model_list = g_slist_reverse(model_list);
for (list=model_list ; list ; list=g_slist_next(list))
  {
  temp = list->data;

/* NEW - process the library symbols */
  for (flist=temp->cores ; flist ; flist=g_slist_next(flist))
    {
    core = flist->data;

    text = g_hash_table_lookup(types, core->atom_label); 

    if (text)
      {
      g_free(core->atom_type);
      core->atom_type = g_strdup(text);
      }
    }

/* init for display (if not already displayed) */
  if (!temp->grafted && !pflag)
    {
    model_prep(temp);
    }
  }

/* CURRENT */
/*
dump_phonon(289, 873, data);
*/

g_hash_table_destroy(types);

fclose(fp);
return(0);
}

/**************************************/
/* show the current file stream error */
/**************************************/
void print_file_error(FILE *fp)
{

printf(">>> [offet = %ld] ", ftell(fp));

if (feof(fp))
  printf("error = EOF\n");
else
  printf("error = %d\n", ferror(fp));

clearerr(fp);
}

/***********************************/
/* GULP dynamics trajectory header */
/***********************************/
/* macro for fortran start/end record padding */
/*
int int_buffer=0;
#define READ_RECORD fread(&int_buffer, sizeof(int), 1, fp)
#define WRITE_RECORD fwrite(&int_buffer, sizeof(int), 1, fp)
*/

/* NB: I experimented with fseek AND fsetpos, but these seem to fail */
/* with EOF errors even tho ftell reports a position < file size */

#define DEBUG_TRJ_HEADER 0
gint read_trj_header(FILE *fp, struct model_pak *model)
{
int num_atoms, periodic;
double version;

/* record 1 - version */
READ_RECORD;
fread(&version, sizeof(version), 1, fp);
if (fabs(version) > 100.0)
  {
  swap_bytes(&version, sizeof(version));
  if (fabs(version) > 100.0)
    {
    printf("Error: file seems garbled.\n");
    return(1);
    }
  model->trj_swap = TRUE;
  }
READ_RECORD;
/* record 2 */
READ_RECORD;
/* # of atoms + shells */
fread(&num_atoms, sizeof(num_atoms), 1, fp);
if (model->trj_swap)
  swap_bytes(&num_atoms, sizeof(num_atoms));
/* dimension */
fread(&periodic, sizeof(periodic), 1, fp);
if (model->trj_swap)
  swap_bytes(&periodic, sizeof(periodic));
READ_RECORD;

#if DEBUG_TRJ_HEADER
printf("trj version: %lf\n", version);
printf(" Swap bytes: %d\n", model->trj_swap);
printf("total atoms: %d\n", num_atoms);
#endif

return(0);
}

/**********************************/
/* GULP dynamics trajectory frame */
/**********************************/
#define DEBUG_TRJ_FRAME 0
gint read_trj_frame(FILE *fp, struct model_pak *model, gint replace_coords)
{
gint i, j;
int num_atoms, periodic;
double time, ke, pe, temp, *x[6];
double cell[9], velc[6];
GSList *list;
struct core_pak *core;
struct shel_pak *shell;

/* standard frame read */
num_atoms = model->expected_cores + model->expected_shells;
periodic = model->periodic;

/* alloc for x,y,z & vx, vy, vz */
for (j=0 ; j<6 ; j++)
  x[j] = g_malloc(num_atoms * sizeof(double));

/* read one frame */
READ_RECORD;

if (!fread(&time, sizeof(time), 1, fp))
  print_file_error(fp);
if (model->trj_swap)
  swap_bytes(&time, sizeof(time));

if (!fread(&ke, sizeof(ke), 1, fp))
  print_file_error(fp);
if (model->trj_swap)
  swap_bytes(&ke, sizeof(ke));

if (!fread(&pe, sizeof(pe), 1, fp))
  print_file_error(fp);
if (model->trj_swap)
  swap_bytes(&pe, sizeof(pe));

if (!fread(&temp, sizeof(temp), 1, fp))
  print_file_error(fp);
if (model->trj_swap)
  swap_bytes(&temp, sizeof(temp));

READ_RECORD;

#if DEBUG_TRJ_FRAME
printf("[%1dD][%d atoms]", periodic, num_atoms);
printf("[t : %.2lf]",time);
printf("[ke : %.2lf]",ke);
printf("[pe : %.2lf]",pe);
printf("[T : %.2lf][%d]\n",temp, replace_coords);
#endif

/* record dynamics frame */
model->gulp.frame_time = time;
model->gulp.frame_ke = ke;
model->gulp.frame_pe = pe;
model->gulp.frame_temp = temp;

/* loop over x,y,z & assoc velocity components */
for (j=0 ; j<6 ; j++)
  {
/* NB: cores first, then shells */
  READ_RECORD;
  for (i=0 ; i<num_atoms ; i++)
    if (!fread(x[j]+i, sizeof(double), 1, fp))
      {
      print_file_error(fp);
      break;
      }
  READ_RECORD;
  }

#if DEBUG_TRJ_FRAME
printf("expected: cores = %d, shells = %d\n",
        model->expected_cores, model->expected_shells);
printf("existing: cores = %d, shells = %d\n",
        g_slist_length(model->cores), g_slist_length(model->shels));
#endif

/* read core data */
list = model->cores;
if (replace_coords)
  {
  for (i=0 ; i<model->expected_cores ; i++)
    {
    if (list)
      core = list->data;
    else
      {
      core = new_core("X", model);
      model->cores = g_slist_append(model->cores, core);
      list = g_slist_last(model->cores);
      }

    core->x[0] = *(x[0]+i);
    core->x[1] = *(x[1]+i);
    core->x[2] = *(x[2]+i);

    core->v[0] = *(x[3]+i);
    core->v[1] = *(x[4]+i);
    core->v[2] = *(x[5]+i);
  
    list = g_slist_next(list);
    }
  }

/* read shells data */
list = model->shels;
if (replace_coords)
  {
  for (i=model->expected_cores ; i<num_atoms ; i++)
    {
    if (list)
      shell = list->data;
    else
      {
      shell = new_shell("X", model);
      model->shels = g_slist_append(model->shels, shell);
      list = g_slist_last(model->shels);
      }

    shell->x[0] = *(x[0]+i);
    shell->x[1] = *(x[1]+i);
    shell->x[2] = *(x[2]+i);
  
    shell->v[0] = *(x[3]+i);
    shell->v[1] = *(x[4]+i);
    shell->v[2] = *(x[5]+i);
  
    list = g_slist_next(list);
    }
  }

/* read cell info */
if (model->gulp.ensemble == NPT)
  {
/* get cell vectors */
  READ_RECORD;
  fread(cell, sizeof(double), 9, fp);
  READ_RECORD;

/* get cell velocities */
  READ_RECORD;
/* calculate nstrains [0, 6] */
  switch (model->periodic)
    {
    case 3:
      j=6;
      break;
    case 2:
      j=3;
      break;
    case 1:
      j=1;
      break;
    default:
      j=0;
    }
  fread(velc, sizeof(double), j, fp);
  READ_RECORD;

/* compute latmat/ilatmat */
  memcpy(model->latmat, cell, 9*sizeof(gdouble));
  model->construct_pbc = TRUE;
  matrix_lattice_init(model); 
  }

/* convert input cartesian to fractional */
if (replace_coords)
  if (model->fractional)
    coords_make_fractional(model);

/* clean up */
for (j=0 ; j<6 ; j++)
  g_free(x[j]);

return(0);
}

/***********************************/
/* GULP dynamics trajectory header */
/***********************************/
gint write_trj_header(FILE *fp, struct model_pak *model)
{
int num_atoms, periodic;
double version = 1.4;

WRITE_RECORD;
fwrite(&version, sizeof(version), 1, fp);
WRITE_RECORD;

WRITE_RECORD;
num_atoms = g_slist_length(model->cores) + g_slist_length(model->shels);
fwrite(&num_atoms, sizeof(num_atoms), 1, fp);
periodic = model->periodic;
fwrite(&periodic, sizeof(periodic), 1, fp);
WRITE_RECORD;

return(0);
}

/****************************************/
/* write GULP dynamics trajectory frame */
/****************************************/
gint write_trj_frame(FILE *fp, struct model_pak *model)
{
gint i, j;
double temp, x[3];
GSList *list;
struct core_pak *core;
struct shel_pak *shell;

/* fake time, ke, pe, and temp */
WRITE_RECORD;
temp = 0.0;
fwrite(&temp, sizeof(double), 4, fp);
WRITE_RECORD;

/* core/shell coordinate write */
for (i=0 ; i<3 ; i++)
  {
  WRITE_RECORD;
  for (list=model->cores ; list ; list=g_slist_next(list))
    {
    core = list->data;

    ARR3SET(x, core->x);
    vecmat(model->latmat, x);
  
    fwrite(&x[i], sizeof(double), 1, fp);
    }
  for (list=model->shels ; list ; list=g_slist_next(list))
    {
    shell = list->data;

    ARR3SET(x, shell->x);
    vecmat(model->latmat, x);
  
    fwrite(&x[i], sizeof(double), 1, fp);
    }
  WRITE_RECORD;
  }

/* core/shell velocity write */
for (i=0 ; i<3 ; i++)
  {
  WRITE_RECORD;
  for (list=model->cores ; list ; list=g_slist_next(list))
    {
    core = list->data;

    x[i] = core->v[i];
  
    fwrite(&x[i], sizeof(double), 1, fp);
    }
  for (list=model->shels ; list ; list=g_slist_next(list))
    {
    shell = list->data;

    x[i] = shell->v[i];
  
    fwrite(&x[i], sizeof(double), 1, fp);
    }
  WRITE_RECORD;
  }

/* lattice vector write */
if (model->gulp.ensemble == NPT)
  {
  WRITE_RECORD;
/* TODO - check matrix storage order is ok (ie in rows or columns?) */
  for (i=0 ; i<9 ; i++)
    {
    temp = model->latmat[i];
    fwrite(&temp, sizeof(double), 1, fp);
    }
  WRITE_RECORD;

/* calculate nstrains [0, 6] */
  switch (model->periodic)
    {
    case 3:
      j=6;
      break;
    case 2:
      j=3;
      break;
    case 1:
      j=1;
      break;
    default:
      j=0;
    }

/* fake lattice velocity write */
  WRITE_RECORD;
  temp = 0.0;
  for (i=0 ; i<j ; i++)
    fwrite(&temp, sizeof(double), 1, fp);
  WRITE_RECORD;
  }

return(0);
}