/*
Copyright (C) 2003 by Andrew Lloyd Rohl
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 <stdio.h>
#include <string.h>
#include "gdis.h"
#include "coords.h"
#include "model.h"
#include "file.h"
#include "parse.h"
#include "interface.h"
/* main structures */
extern struct sysenv_pak sysenv;
extern struct elem_pak elements[];
/* globals */
gint optcell;
/*******************************************/
/* read single ABINIT output configuration */
/*******************************************/
#define DEBUG_READ_ABOUT 0
gint read_about_block(FILE *fp, struct model_pak *model)
{
gint i, num_tokens;
gchar **buff, line[LINELEN];
gdouble acell[3];
GString *title;
GSList *clist;
struct core_pak *core;
clist = model->cores;
if (optcell > 0)
{
/* go to acell line and read it */
if (fgetline(fp, line))
{
gui_text_show(ERROR, "unexpected end of file reading cell dimensions in animation\n");
return(2);
}
buff = tokenize(line, &num_tokens);
acell[0] = str_to_float(*(buff+1)) * AU2ANG;
acell[1] = str_to_float(*(buff+2)) * AU2ANG;
acell[2] = str_to_float(*(buff+3)) * AU2ANG;
g_strfreev(buff);
/* get the cell vectors i.e. rprim.acell */
/* NB: gdis wants transposed ABINIT matrix */
for (i=0; i<3; i++)
{
if (fgetline(fp, line))
{
gui_text_show(ERROR, "unexpected end of file reading cell dimensions\n");
return(2);
}
buff = tokenize(line + 8, &num_tokens);
model->latmat[0+i] = acell[i] * str_to_float(*(buff+0));
model->latmat[3+i] = acell[i] * str_to_float(*(buff+1));
model->latmat[6+i] = acell[i] * str_to_float(*(buff+2));
g_strfreev(buff);
}
for (i=0; i<4; i++)
{
if (fgetline(fp, line))
{
gui_text_show(ERROR, "unexpected end of file skipping to coordinates in animation\n");
return(2);
}
}
}
/* get 1st line of coords */
if (fgetline(fp, line))
return(1);
buff = tokenize(line, &num_tokens);
while (g_ascii_strncasecmp(line, " Cartesian forces", 17) != 0)
{
if (clist)
{
core = (struct core_pak *) clist->data;
clist = g_slist_next(clist);
}
else
{
core = new_core(*(buff+4), model);
model->cores = g_slist_append(model->cores, core);
}
core->x[0] = AU2ANG*str_to_float(*(buff+0));
core->x[1] = AU2ANG*str_to_float(*(buff+1));
core->x[2] = AU2ANG*str_to_float(*(buff+2));
#if DEBUG_READ_ABOUT
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))
return(2);
buff = tokenize(line, &num_tokens);
}
/* get gradients */
model->abinit.max_grad = str_to_float(*(buff+4));
model->abinit.rms_grad = str_to_float(*(buff+5));
g_strfreev(buff);
/* get energy */
if (fgetline(fp, line))
return(2);
while (g_ascii_strncasecmp(line, " At the end of Broyden", 22) != 0)
{
if (fgetline(fp, line))
return(2);
}
buff = tokenize(line, &num_tokens);
model->abinit.energy = str_to_float(*(buff+9));
title = g_string_new("");
g_string_append_printf(title, "E");
g_string_append_printf(title, " = %.4f Ha, ", model->abinit.energy);
g_string_append_printf(title, "max grad = %.5f", model->abinit.max_grad);
model->title = g_strdup(title->str);
g_string_free(title, TRUE);
return(0);
}
/*******************************/
/* ABINIT output frame reading */
/*******************************/
gint read_about_frame(FILE *fp, struct model_pak *model)
{
return(read_about_block(fp, model));
}
/********************************/
/* Read in a ABINIT output file */
/********************************/
gint read_about(gchar *filename, struct model_pak *model)
{
gint i, flag, frame, num_tokens, tot_tokens;
gint natom=0, ntype=0, *types;
GString *title;
gchar **buff, **curr_token, *ptr, line[LINELEN];
struct core_pak *core;
GSList *clist;
FILE *fp;
fp = fopen(filename, "rt");
if (!fp)
return(1);
flag=frame=0;
while (!fgetline(fp, line))
{
/* space group info */
if (g_ascii_strncasecmp(" Symmetries :", line, 13) == 0)
{
ptr = g_strrstr(line, "(#");
model->sginfo.spacenum = (gint) str_to_float(ptr+2);
}
/* default cell dimensions */
/* NB: gdis wants transposed ABINIT matrix */
if (g_ascii_strncasecmp(" Real(R)+Recip(G)", line, 17) == 0)
{
for (i=0; i<3; i++)
{
if (fgetline(fp, line))
{
gui_text_show(ERROR, "unexpected end of file reading cell dimensions\n");
return(2);
}
buff = tokenize(line, &num_tokens);
model->latmat[0+i] = str_to_float(*(buff+1))*AU2ANG;
model->latmat[3+i] = str_to_float(*(buff+2))*AU2ANG;
model->latmat[6+i] = str_to_float(*(buff+3))*AU2ANG;
g_strfreev(buff);
}
model->construct_pbc = TRUE;
model->periodic = 3;
/* last part of data in output file before minimisation */
break;
}
/* optimisation type */
if (g_ascii_strncasecmp(" optcell", line, 10) == 0)
{
buff = tokenize(line, &num_tokens);
optcell = (gint) str_to_float(*(buff+1));
g_strfreev(buff);
}
/* coordinates */
if (g_ascii_strncasecmp(" natom", line, 10) == 0)
{
buff = tokenize(line, &num_tokens);
natom = (gint) str_to_float(*(buff+1));
g_strfreev(buff);
}
if (g_ascii_strncasecmp(" ntype", line, 10) == 0)
{
buff = tokenize(line, &num_tokens);
ntype = (gint) str_to_float(*(buff+1));
g_strfreev(buff);
}
if (g_ascii_strncasecmp(" type", line, 10) == 0)
{
tot_tokens = 0;
buff = tokenize(line, &num_tokens);
curr_token = buff + 1;
while (tot_tokens < natom)
{
if (*curr_token == NULL)
{
g_strfreev(buff);
buff = get_tokenized_line(fp, &num_tokens);
curr_token = buff;
}
core = new_core("X", model);
model->cores = g_slist_prepend(model->cores, core);
core->atom_code = (gint) str_to_float(*curr_token) - 1;
tot_tokens++;
curr_token++;
}
model->cores = g_slist_reverse(model->cores);
g_strfreev(buff);
flag++;
}
if (g_ascii_strncasecmp(" xangst", line, 10) == 0)
{
clist = model->cores;
if (clist == NULL)
{
gui_text_show(ERROR, "no atoms found\n");
return(2);
}
buff = tokenize(line+10, &num_tokens);
for (i=0; i<natom; i++)
{
core = (struct core_pak *) clist->data;
core->x[0] = str_to_float(*(buff+0));
core->x[1] = str_to_float(*(buff+1));
core->x[2] = str_to_float(*(buff+2));
g_strfreev(buff);
buff = get_tokenized_line(fp, &num_tokens);
clist = g_slist_next(clist);
}
}
if (g_ascii_strncasecmp(" znucl", line, 10) == 0)
{
types = g_malloc(ntype * sizeof(gint));
i = 0;
tot_tokens = 0;
buff = tokenize(line, &num_tokens);
curr_token = buff + 1;
while (tot_tokens < ntype)
{
if (*curr_token == NULL)
{
g_strfreev(buff);
buff = get_tokenized_line(fp, &num_tokens);
curr_token = buff;
}
types[i] = (gint) str_to_float(*curr_token);
i++;
tot_tokens++;
curr_token++;
}
g_strfreev(buff);
/* now update cores to get right element type */
clist = model->cores;
if (clist == NULL)
{
gui_text_show(ERROR, "no atoms found\n");
return(2);
}
for (i=0; i<natom; i++)
{
core = clist->data;
core_init(elements[types[core->atom_code]].symbol, core, model);
clist = g_slist_next(clist);
}
}
}
/* Additional info */
while (!fgetline(fp, line))
{
/* energy */
if (g_ascii_strncasecmp(line, " Ewald =", 9) == 0)
{
printf("reading energy\n");
model->abinit.energy = str_to_float(&line[51]);
printf("energy is %lf\n", model->abinit.energy);
}
/* gradient */
else if ((g_ascii_strncasecmp(line, " frms,max,avg=", 14) == 0) && g_strrstr(line, "h/b"))
{
printf("reading gradient\n");
buff = tokenize(line, &num_tokens);
model->abinit.max_grad = str_to_float(*(buff+2));
model->abinit.rms_grad = str_to_float(*(buff+1));
g_strfreev(buff);
printf("gradient is %lf\n", model->abinit.max_grad);
}
/* coordinates */
else if ((g_ascii_strncasecmp(line, " Cartesian coordinates (bohr)", 29) == 0 && optcell == 0) || (g_ascii_strncasecmp(line, " Unit cell characteristics :", 28) == 0 && optcell > 0))
{
/* go through all frames to count them */
add_frame_offset(fp, model);
read_about_block(fp, model);
frame++;
}
}
/* done */
if (flag)
{
strcpy(model->filename, filename);
g_free(model->basename);
model->basename = strdup_basename(filename);
model->num_frames = model->cur_frame = frame;
model->cur_frame--;
title = g_string_new("");
g_string_append_printf(title, "E");
g_string_append_printf(title, " = %.4f Ha, ", model->abinit.energy);
g_string_append_printf(title, "max grad = %.5f", model->abinit.max_grad);
model->title = g_strdup(title->str);
g_string_free(title, TRUE);
model_prep(model);
}
else
return(2);
return(0);
}
syntax highlighted by Code2HTML, v. 0.9.1