/*
Copyright (C) 2004 by Andrew Lloyd Rohl

andrew@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 "gdis.h"
#include "coords.h"
#include "file.h"
#include "parse.h"
#include "model.h"
#include "interface.h"

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

/*******************************************/
/* read single CASTEP output configuration */
/*******************************************/
#define DEBUG_READ_CASTEP_OUT 0
gint read_castep_out_block(FILE *fp, struct model_pak *model)
{
gint i;
gint num_tokens;
gint no_grad = 0;
gint last_frame = FALSE;
gdouble grad, max_grad = 0.0, rms_grad = 0.0;
gdouble pressure;
gchar **buff, line[LINELEN], *ext, *text;
GString *title, *grad_string, *pressure_string;
GSList *clist;
struct core_pak *core;

clist = model->cores;

/* read to end of iteration */
while (TRUE)
  {
  if (fgetline(fp, line))
    {
    gui_text_show(ERROR, "unexpected end of file reading to end of iteration\n");
    return(2);
    }
    
  if (g_ascii_strncasecmp(line, "==========", 10) == 0)
      break;
      
  if (g_ascii_strncasecmp(line, "Writing model to", 16) == 0)
      break;
  
/* read cell vectors? */
  if (g_strrstr(line, "Real Lattice") != NULL)
    {
    /* read in cell vectors */
    /* NB: gdis wants transposed matrix */
    for (i=0; i<3; i++)
      {
      if (fgetline(fp, line))
        {
        gui_text_show(ERROR, "unexpected end of file reading cell vectors\n");
        return(2);
        }
      buff = tokenize(line, &num_tokens);
      model->latmat[0+i] = str_to_float(*(buff+0));
      model->latmat[3+i] = str_to_float(*(buff+1));
      model->latmat[6+i] = str_to_float(*(buff+2));
      g_strfreev(buff);
      }
    }
 
/* read coordinates */
  if (g_strrstr(line, "x-") != NULL)
    {
    if (fgetline(fp, line))
      {
      gui_text_show(ERROR, "unexpected end of file skipping to coordinates\n");
      return(2);
      }

    buff = tokenize(line, &num_tokens);
    while (num_tokens == 7)
      {
      if (clist)
        {
        core = (struct core_pak *) clist->data;
        clist = g_slist_next(clist);
        }
      else
        {
        core = new_core(*(buff+1), model);
        model->cores = g_slist_append(model->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));
      #if DEBUG_READ_CASTEP_OUT
      printf("new coords %f %f %f\n", core->x[0],  core->x[1], core->x[2]);
      #endif

    /* get next line */
      g_strfreev(buff);
      if (fgetline(fp, line))
        {
        gui_text_show(ERROR, "unexpected end of file reading coordinates\n");
        return(2);
        }
      buff = tokenize(line, &num_tokens);
      }
    g_strfreev(buff);
    clist = model->cores;
    }

/* SCF or Final Energy */
  if (g_strrstr(line, "Number of kpoints used") != NULL)
    {
    buff = g_strsplit(line, "=", 2);
    g_strchug(g_strchomp(*(buff+1)));
    property_add_ranked(10, "K-points", *(buff+1), model);
    g_strfreev(buff);
    }

  if (g_strrstr(line, "Files used for pseudopotentials:") != NULL)
    {
    if (fgetline(fp, line))
      {
      gui_text_show(ERROR, "unexpected end of file reading pseudotential file names\n");
      return(2);
      }
      buff = tokenize(line, &num_tokens);
      ext = find_char(*(buff+1), '.', LAST);
      ext++;
      if (g_strrstr(ext, "usp") != NULL)
        property_add_ranked(10, "Pseudopotentials", "Ultrasoft", model);
      else if (g_strrstr(ext, "recpot") != NULL)
        property_add_ranked(10, "Pseudopotentials", "Norm conserving", model);
      else
        property_add_ranked(10, "Pseudopotentials", "Unknown", model);
    }
  if (g_ascii_strncasecmp(line, "Final energy", 12) == 0)
    {
    buff = g_strsplit(line, "=", 2);
    text = format_value_and_units(*(buff+1), 5);
    property_add_ranked(3, "Energy", text, model);
    g_free(text);
    model->castep.energy = str_to_float(*(buff+1));
    model->castep.have_energy = TRUE;
    g_strfreev(buff);
    }
  if (g_strrstr(line, "Final Enthalpy") != NULL)
    {
    last_frame = TRUE;
    model->castep.min_ok = TRUE;
    buff = g_strsplit(line, "=", 2);
    text = format_value_and_units(*(buff+1), 5);
    property_add_ranked(3, "Energy", text, model);
    g_free(text);
    model->castep.energy = str_to_float(*(buff+1));
    model->castep.have_energy = TRUE;
    g_strfreev(buff);
    }

  /* gradients */
  if (g_strrstr(line, "* -") != NULL)
    {
    for (i=0; i<2; i++)
      {
      if (fgetline(fp, line))
        {
        gui_text_show(ERROR, "unexpected end of file skipping to gradients\n");
        return(2);
        }
      }
    /* read gradients */
    if (fgetline(fp, line))
      {
      gui_text_show(ERROR, "unexpected end of file reading gradients\n");
      return(2);
      }
    buff = tokenize(line, &num_tokens);
    while (num_tokens == 7)
      {
      for (i=0; i<3; i++)
        {
        grad = str_to_float(*(buff+3+i));
        no_grad++;
        rms_grad = grad*grad;
        max_grad = MAX(fabs(grad), max_grad);
      }
      /* get next line */
      g_strfreev(buff);
      if (fgetline(fp, line))
        {
        gui_text_show(ERROR, "unexpected end of file reading gradients\n");
        return(2);
        }
      buff = tokenize(line, &num_tokens);
      }
    g_strfreev(buff);
    model->castep.rms_grad = sqrt(rms_grad/no_grad);
    model->castep.max_grad = max_grad;
    grad_string = g_string_new("");
    g_string_append_printf(grad_string, "%.5f %s", model->castep.rms_grad, property_lookup("Gradient Units", model));
    property_add_ranked(4, "RMS Gradient", grad_string->str, model);
    g_string_free(grad_string, TRUE); 

    #if DEBUG_READ_CASTEP_OUT
    printf("RMS gradient: %f Max gradient: %f\n", model->castep.rms_grad, model->castep.max_grad);
    #endif
    }
      
  if (g_strrstr(line, "Pressure") != NULL)
    {
    /* read pressure */
    buff = tokenize(line, &num_tokens);
    pressure = str_to_float(*(buff+2));
    pressure_string = g_string_new("");
    g_string_append_printf(pressure_string, "%.4f %s", pressure, property_lookup("Pressure Units", model));
    property_add_ranked(5, "Pressure", pressure_string->str, model);
    g_string_free(pressure_string, TRUE);
    g_strfreev(buff);
    }
  }
  
g_free(model->title);
title = g_string_new("");
g_string_append_printf(title, "E");
g_string_append_printf(title, "(DFT)");
/* TODO read energy units from CASTEP file */
g_string_append_printf(title, " = %.5f eV", model->castep.energy);
g_string_append_printf(title, ", grad = %.5f", model->castep.rms_grad);
model->title = g_strdup(title->str);
g_string_free(title, TRUE); 

return(0);
}

/*******************************/
/* CASTEP output frame reading */
/*******************************/
gint read_castep_out_frame(FILE *fp, struct model_pak *model)
{
    return(read_castep_out_block(fp, model));
}

/********************************/
/* Read in a CASTEP output file */
/********************************/
gint read_castep_out(gchar *filename, struct model_pak *model)
{
gint frame;
gchar **buff, line[LINELEN], *text;
FILE *fp;

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

model->construct_pbc = TRUE;
model->periodic = 3;
model->fractional = TRUE;
  
frame=0;
while (!fgetline(fp, line))
  {
  /* find relevent units */
  if (g_strrstr(line, "force unit") != NULL)
    {
    buff = g_strsplit(line, ":", 2);
    g_strstrip(*(buff+1));
    property_add_ranked(0, "Gradient Units", *(buff+1), model);
    g_strfreev(buff);
    }
  if (g_strrstr(line, "pressure unit") != NULL)
    {
    buff = g_strsplit(line, ":", 2);
    g_strstrip(*(buff+1));
    property_add_ranked(0, "Pressure Units", *(buff+1), model);
    g_strfreev(buff);
    }
  /* store selected parameters in hash table */
  if (g_ascii_strncasecmp(line, " type of calculation", 20) == 0)
    {
    buff = g_strsplit(line, ":", 2);
    g_strstrip(*(buff+1));
    if (g_ascii_strncasecmp(*(buff+1), "geometry optimization", 21) == 0)
      property_add_ranked(2, "Calculation", "Optimisation", model);
    else if (g_ascii_strncasecmp(*(buff+1), "single point energy", 19) == 0)
      property_add_ranked(2, "Calculation", "Single Point", model);
    else
      property_add_ranked(2, "Calculation", *(buff+1), model);
    g_strfreev(buff);
    }
  if (g_ascii_strncasecmp(line, " using functional", 17) == 0)
    {
    buff = g_strsplit(line, ":", 2);
    g_strstrip(*(buff+1));
    if (g_ascii_strncasecmp(*(buff+1), "Perdew Burke Ernzerhof", 22) == 0)
      property_add_ranked(7, "Functional", "PBE", model);
    else if (g_ascii_strncasecmp(*(buff+1), "Perdew Wang (1991)", 18) == 0)
      property_add_ranked(7, "Functional", "PW91", model);
    else
      property_add_ranked(7, "Functional", *(buff+1), model);
    g_strfreev(buff);
    }
  if (g_ascii_strncasecmp(line, " plane wave basis set cut-off", 29) == 0)
    {
    buff = g_strsplit(line, ":", 2);
    g_strstrip(*(buff+1));
    text = format_value_and_units(*(buff+1), 1);
    property_add_ranked(6, "Basis Cutoff", text, model);
    g_free(text);
    }  
  if (g_ascii_strncasecmp(line, " size of standard grid", 22) == 0)
    {
    buff = g_strsplit(line, ":", 2);
    g_strstrip(*(buff+1));
    text = format_value_and_units(*(buff+1), 2);
    property_add_ranked(10, "Grid", text, model);
    g_free(text);
    }  
 
  /* read coordinates */
	if ((g_strrstr(line, "Unit Cell") != NULL) || (g_strrstr(line, "Cell Contents") != NULL))
		{
		/* go through all frames to count them */
		add_frame_offset(fp, model);
		read_castep_out_block(fp, model);
		frame++;
    }
  }

/* done */
strcpy(model->filename, filename);
g_free(model->basename);
model->basename = strdup_basename(filename);
model->num_frames = model->cur_frame = frame;
model->cur_frame--;

model_prep(model);

return(0);
}