/*
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 <stdio.h>
#include <string.h>
#include <time.h>
#include "gdis.h"
#include "coords.h"
#include "file.h"
#include "parse.h"
#include "matrix.h"
#include "model.h"
#include "interface.h"
#define DEBUG_MORE 0
#define MAX_KEYS 15
/* main structures */
extern struct sysenv_pak sysenv;
extern struct elem_pak elements[];
/***************/
/* CIF writing */
/***************/
gint write_cif(gchar *filename, struct model_pak *data)
{
gint flag=0;
gdouble tmat[9], x[3], depth=1.0;
GSList *list;
FILE *fp;
time_t t1;
struct core_pak *core;
/* init & checks */
g_return_val_if_fail(data != NULL, 1);
fp = fopen(filename, "wt");
g_return_val_if_fail(fp != NULL, 2);
/* is it a surface with negative z? */
if (data->periodic == 2)
{
if (g_slist_length(data->cores))
{
core = (struct core_pak *) g_slist_nth_data(data->cores, 0);
if (core->x[2] < 0.0)
flag++;
}
}
/* rot to make z positive */
matrix_rotation(&tmat[0], PI, PITCH);
t1 = time(NULL);
fprintf(fp, "data_block_1\n");
fprintf(fp, "_audit_creation_date '%s'\n", g_strstrip(ctime(&t1)));
fprintf(fp, "_audit_creation_method 'generated by GDIS v%4.2f'\n", VERSION);
fprintf(fp, "\n\n");
if (data->periodic)
{
fprintf(fp, "_cell_length_a %8.4f\n",data->pbc[0]);
fprintf(fp, "_cell_length_b %8.4f\n",data->pbc[1]);
if (data->periodic == 2)
{
/* get depth info */
depth = (data->surface.region[0]+data->surface.region[1])*data->surface.depth;
/* no depth info - make it large enough to fit everything */
if (depth < POSITION_TOLERANCE)
depth = 2.0*data->rmax;
fprintf(fp, "_cell_length_c %8.4f\n", depth);
}
else
fprintf(fp, "_cell_length_c %8.4f\n",data->pbc[2]);
fprintf(fp, "_cell_angle_alpha %8.2f\n",180.0*data->pbc[3]/PI);
fprintf(fp, "_cell_angle_beta %8.2f\n",180.0*data->pbc[4]/PI);
fprintf(fp, "_cell_angle_gamma %8.2f\n",180.0*data->pbc[5]/PI);
fprintf(fp, "\n\n");
fprintf(fp, "_symmetry_space_group_name_H-M '%s'\n",data->sginfo.spacename);
fprintf(fp, "_symmetry_Int_Tables_number %d\n",data->sginfo.spacenum);
fprintf(fp, "\n\n");
}
/* coords - format section */
fprintf(fp, "loop_\n");
fprintf(fp, "_atom_site_label\n");
if (data->periodic)
{
fprintf(fp, "_atom_site_fract_x\n");
fprintf(fp, "_atom_site_fract_y\n");
fprintf(fp, "_atom_site_fract_z\n");
}
else
{
fprintf(fp, "_atom_site_cartn_x\n");
fprintf(fp, "_atom_site_cartn_y\n");
fprintf(fp, "_atom_site_cartn_z\n");
}
/* coords - data section */
for (list=data->cores ; list ; list=g_slist_next(list))
{
core = list->data;
/* only the asymmetric cell if periodic */
if (data->periodic && !core->primary)
continue;
/* NB: want fractional if 3D periodic, otherwise cartesian */
ARR3SET(x, core->x);
/* transformation needed? */
if (flag)
vecmat(tmat, x);
/* make the z coordinate "fractional" for a surface */
if (data->periodic == 2)
x[2] /= depth;
fprintf(fp, "%2s %9.4f %9.4f %9.4f\n",
elements[core->atom_code].symbol, x[0], x[1], x[2]);
}
fprintf(fp, "\n\n");
fclose(fp);
return(0);
}
/***************/
/* CIF loading */
/***************/
#define DEBUG_LOAD_CIF 0
gint read_cif(gchar *filename, struct model_pak *data)
{
gint i, j, n, first, new, order, pos, keyword, loop_count=0;
gint min_tokens, num_tokens, len, flag;
gint s_pos, l_pos, x_pos, y_pos, z_pos, o_pos;
gchar **buff, *tmp, *name=NULL, *line;
gdouble sign, mat[9], off[3];
GSList *list=NULL;
struct core_pak *core;
FILE *fp;
/* checks */
g_return_val_if_fail(data != NULL, 1);
g_return_val_if_fail(filename != NULL, 2);
fp = fopen(filename, "rt");
if (!fp)
return(3);
/* atom counter */
n=-1;
new = first = data->number;
data->id = -1;
for(;;)
{
/* end if EOF */
line = file_read_line(fp);
if (!line)
break;
/*
if (fgetline(fp, line))
break;
*/
/* search for data */
list = get_keywords(line);
cif_parse:;
if (list != NULL)
{
keyword = GPOINTER_TO_INT(list->data);
switch(keyword)
{
/* model labels */
/* FIXME - needs to search for the 1st occurrance of ' or " & then get the string */
/*
case CIF_CHEMICAL_NAME:
tmp = g_strdup(get_token_pos(line,1));
for (i=0 ; i<strlen(tmp) ; i++)
if (*(tmp+i) == '\'')
*(tmp+i) = ' ';
g_free(data->basename);
data->basename = g_strdup(g_strstrip(tmp));
g_free(tmp);
break;
case CIF_MINERAL_NAME:
tmp = g_strdup(get_token_pos(line,1));
for (i=0 ; i<strlen(tmp) ; i++)
if (*(tmp+i) == '\'')
*(tmp+i) = ' ';
g_free(data->basename);
data->basename = g_strdup(g_strstrip(tmp));
g_free(tmp);
break;
*/
/* stopgap model name */
case CIF_DATA_START:
name = g_strdup(g_strstrip(&line[5]));
break;
/* candidate new model trigger */
/* CIF is a pain - there seems to be no clear cut new structure */
/* trigger (in a multi structure file) that everyone uses */
case CIF_AUDIT:
/* skip this the 1st time (we've already alloc'd a data pointer for safety) */
new++;
if ((new-first) > 1)
{
/* NEW - some dodgy models have no atom data - so allow new model */
/* creation if we have (the bare minimum) some new cell parameters */
if (!data->periodic)
{
/* no cell parameters found yet - don't create a new model */
new--;
break;
}
# if DEBUG_LOAD_CIF
printf("Found %d atoms. [reset]\n", n);
#endif
/* alloc new pointer */
data = model_new();
if (data == NULL)
goto cif_done;
}
/* if we found a name - assign it now */
if (name)
{
g_free(data->basename);
data->basename = name;
name = NULL;
}
#if DEBUG_LOAD_CIF
printf("Start of new model %d: %s.\n",data->number,data->basename);
#endif
break;
/* model specific data - the *data ptr MUST be allocated */
case CIF_CELL_A:
buff = get_tokens(line, 3);
data->pbc[0] = str_to_float(*(buff+1));
data->periodic++;
g_strfreev(buff);
break;
case CIF_CELL_B:
buff = get_tokens(line, 3);
data->pbc[1] = str_to_float(*(buff+1));
data->periodic++;
g_strfreev(buff);
break;
case CIF_CELL_C:
buff = get_tokens(line, 3);
data->pbc[2] = str_to_float(*(buff+1));
data->periodic++;
g_strfreev(buff);
break;
case CIF_CELL_ALPHA:
buff = get_tokens(line, 3);
data->pbc[3] = PI*str_to_float(*(buff+1))/180.0;
g_strfreev(buff);
break;
case CIF_CELL_BETA:
buff = get_tokens(line, 3);
data->pbc[4] = PI*str_to_float(*(buff+1))/180.0;
g_strfreev(buff);
break;
case CIF_CELL_GAMMA:
buff = get_tokens(line, 3);
data->pbc[5] = PI*str_to_float(*(buff+1))/180.0;
g_strfreev(buff);
break;
case CIF_SPACE_NAME:
/* remove the enclosing ' characters */
tmp = g_strdup(get_token_pos(line,1));
for (i=0 ; i<strlen(tmp) ; i++)
if (*(tmp+i) == '\'')
*(tmp+i) = ' ';
/* store the name, stripping spaces */
data->sginfo.spacename = g_strdup(g_strstrip(tmp));
/* indicate that name should used in lookup */
data->sginfo.spacenum = -1;
#if DEBUG_LOAD_CIF
printf("spacegroup: [%s]\n",data->sginfo.spacename);
#endif
g_free(tmp);
break;
case CIF_SPACE_NUM:
buff = get_tokens(line, 3);
g_strfreev(buff);
break;
case CIF_EQUIV_SITE:
loop_count++;
break;
/* NEW - reinserted the lost symmetry matrix code */
case CIF_EQUIV_POS:
/* allocate for order number of pointers (to matrices) */
data->sginfo.matrix = (gdouble **) g_malloc(sizeof(gdouble *));
data->sginfo.offset = (gdouble **) g_malloc(sizeof(gdouble *));
data->sginfo.order=order=0;
for(;;)
{
/* terminate on EOF, */
g_free(line);
line = file_read_line(fp);
if (!line)
break;
/* blank line, */
g_strstrip(line);
if (!strlen(line))
break;
/* or a new command */
list = get_keywords(line);
if (list)
goto cif_parse;
/* TODO - make this parsing a subroutine */
for(i=0 ; i<strlen(line) ; i++)
if (*(line+i) == '\'' || *(line+i) == ',')
*(line+i) = ' ';
g_strstrip(line);
buff = tokenize(line, &num_tokens);
n = loop_count;
while (n < num_tokens-2)
{
/* FIXME - yet another mess that a linked list would greatly simplify */
/* number of ops */
data->sginfo.matrix = (gdouble **) g_renew
(gdouble *, data->sginfo.matrix , (order+1));
data->sginfo.offset = (gdouble **) g_renew
(gdouble *, data->sginfo.offset , (order+1));
/* actual op */
*(data->sginfo.matrix+order) = (gdouble *) g_malloc(9*sizeof(gdouble));
*(data->sginfo.offset+order) = (gdouble *) g_malloc(3*sizeof(gdouble));
#if DEBUG_LOAD_CIF
printf("[%s] [%s] [%s]\n", *(buff+n+0), *(buff+n+1), *(buff+n+2));
#endif
VEC3SET(&mat[0], 0.0, 0.0, 0.0);
VEC3SET(&mat[3], 0.0, 0.0, 0.0);
VEC3SET(&mat[6], 0.0, 0.0, 0.0);
VEC3SET(&off[0], 0.0, 0.0, 0.0);
for (i=0 ; i<3 ; i++)
{
pos=0;
sign=1.0;
for (j=0 ; j<strlen(*(buff+i+n)) ; j++)
{
switch (*(*(buff+i+n)+j))
{
case '-':
sign = -1.0;
break;
case '+':
sign = +1.0;
break;
case 'x':
mat[i*3] = sign*1.0;
pos++;
break;
case 'y':
mat[i*3 + 1] = sign*1.0;
pos++;
break;
case 'z':
mat[i*3 + 2] = sign*1.0;
pos++;
break;
/* FIXME - a bit crude */
case '/':
g_free(line);
line = g_strndup(*(buff+i+n)+j-1, 3);
off[i] = sign * str_to_float(line);
break;
}
}
}
ARR3SET((*(data->sginfo.matrix+order)+0), &mat[0]);
ARR3SET((*(data->sginfo.matrix+order)+3), &mat[3]);
ARR3SET((*(data->sginfo.matrix+order)+6), &mat[6]);
ARR3SET((*(data->sginfo.offset+order)+0), &off[0]);
#if DEBUG_LOAD_CIF
P3MAT("output: ", *(data->sginfo.matrix+order));
P3VEC("output: ", *(data->sginfo.offset+order));
printf("\n\n");
#endif
order++;
data->sginfo.order++;
n += 3;
}
g_strfreev(buff);
}
#if DEBUG_LOAD_CIF
printf("Found %d symmetry matrices.\n", order);
#endif
data->sginfo.order = order;
break;
case CIF_LOOP_START:
loop_count = 0;
break;
/* parsing for column info ie x,y,z positions frac/cart etc. */
case CIF_ATOM_SITE:
s_pos = l_pos = x_pos = y_pos = z_pos = o_pos = -1;
while (g_strrstr(line, "atom_site") != NULL)
{
if (g_strrstr(line, "type_symbol"))
s_pos = loop_count;
if (g_strrstr(line, "label"))
l_pos = loop_count;
if (g_strrstr(line, "cart_x"))
{
x_pos = loop_count;
data->fractional = FALSE;
}
if (g_strrstr(line, "cart_y"))
{
y_pos = loop_count;
data->fractional = FALSE;
}
if (g_strrstr(line, "cart_z"))
{
z_pos = loop_count;
data->fractional = FALSE;
}
if (g_strrstr(line, "fract_x"))
{
x_pos = loop_count;
data->fractional = TRUE;
}
if (g_strrstr(line, "fract_y"))
{
y_pos = loop_count;
data->fractional = TRUE;
}
if (g_strrstr(line, "fract_z"))
{
z_pos = loop_count;
data->fractional = TRUE;
}
if (g_strrstr(line, "occupancy"))
o_pos = loop_count;
/* get next line and keyword list */
loop_count++;
g_free(line);
line = file_read_line(fp);
if (!line)
goto cif_done;
}
/* either symbol or label can be present & used for atom id purposes */
if (s_pos < 0)
s_pos = l_pos;
if (l_pos < 0)
l_pos = s_pos;
/* check for minimum data */
if (s_pos < 0 || x_pos < 0 || y_pos < 0 || z_pos < 0)
{
#if DEBUG_LOAD_CIF
printf("read_cif() warning: incomplete cif file? [%d:%d:%d:%d]\n",
s_pos, x_pos, y_pos, z_pos);
#endif
break;
}
/* the expected number of tokens */
min_tokens = loop_count;
#if DEBUG_LOAD_CIF
printf(" min tokens: %d\n", min_tokens);
printf("data format: [%d] (%d) - [%d] [%d] [%d] (%d)",
s_pos, l_pos, x_pos, y_pos, z_pos, o_pos);
if (data->fractional)
printf("(frac)\n");
else
printf("(cart)\n");
#endif
/* while no new keywords found */
n = 0;
list = get_keywords(line);
while (list == NULL)
{
buff = tokenize(line, &num_tokens);
/* NB: cif is stupid - it allows data to continue on the next line, */
/* until it gets its number of tokens */
/* hopefully, this'll allow us to get something, even on short lines */
if (num_tokens >= z_pos)
{
/* NEW - ignore labelled (*) symmetry equiv. atoms */
flag = TRUE;
if (l_pos >= 0 && l_pos < num_tokens)
{
len = strlen(*(buff+l_pos));
if ((*(buff+l_pos))[len-1] == '*')
flag = FALSE;
}
if (elem_symbol_test(*(buff+s_pos)) && flag)
{
if (l_pos >= 0 && l_pos < num_tokens)
core = core_new(*(buff+s_pos), *(buff+l_pos), data);
else
core = core_new(*(buff+s_pos), NULL, data);
data->cores = g_slist_prepend(data->cores, core);
core->x[0] = str_to_float(*(buff+x_pos));
core->x[1] = str_to_float(*(buff+y_pos));
core->x[2] = str_to_float(*(buff+z_pos));
/* only get occupancy if we're sure we have enough tokens */
if (o_pos > -1 && num_tokens >= min_tokens)
{
core->sof = str_to_float(*(buff+o_pos));
core->has_sof = TRUE;
}
n++;
}
}
#if DEBUG_LOAD_CIF
else
printf("not enough tokens found.\n");
#endif
/* get next line */
g_strfreev(buff);
g_free(line);
line = file_read_line(fp);
if (!line)
goto cif_done;
list = get_keywords(line);
/* CURRENT - we really want this keyword list parsed again, just in case */
/* a new model trigger (ie audit_creation_date) was found */
if (list)
goto cif_parse;
}
break;
}
}
}
cif_done:;
g_free(line);
/* yet another hack to support the dodgy CIF standard */
j = new - first;
if (!j && n)
{
/* no new model was triggered - but we found atoms, so we'll */
/* assume only one model was in the file & hope for the best */
new++;
}
#if DEBUG_LOAD_CIF
printf("Found %d atoms.\n", n);
printf("Found %d symmetry matrices.\n", data->sginfo.order);
printf("Found %d model(s)\n", new-first);
#endif
/* setup for display */
for (list=sysenv.mal ; list ; list=g_slist_next(list))
{
data = list->data;
if (data->id == -1)
{
data->id = CIF;
data->cores = g_slist_reverse(data->cores);
model_prep(data);
}
}
#if DEBUG_LOAD_CIF
printf("Setup %d model(s).\n",new-first);
#endif
/* clean up & exit */
if (list)
g_slist_free(list);
return(0);
}
syntax highlighted by Code2HTML, v. 0.9.1