/*
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);
}
syntax highlighted by Code2HTML, v. 0.9.1