/*
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 "gdis.h"
#include "coords.h"
#include "file.h"
#include "parse.h"
#include "matrix.h"
#include "model.h"
#include "space.h"
#include "scan.h"
#include "interface.h"
#define DEBUG_MORE 0
#define MAX_KEYS 15
/* main structures */
extern struct sysenv_pak sysenv;
extern struct elem_pak elements[];
/*******************************/
/* MARVIN specific file saving */
/*******************************/
gint write_marvin(gchar *filename, struct model_pak *model)
{
gint mol, region;
gdouble vec[3];
GSList *list;
struct core_pak *core;
struct shel_pak *shel;
FILE *fp;
/* open the file */
fp = fopen(filename,"w");
if (!fp)
return(1);
gui_text_show(STANDARD, "Saving file in MARVIN restart format!\n");
/* save via the number of molecules in the */
/* model - since this allows mol numbering */
mol=0;
/* four region types printing */
for (region=REGION1A ; region<REGION2B ; region++)
{
if (model->region_empty[region])
continue;
switch(region)
{
case REGION1A:
fprintf(fp,"coordinates 1 A\n");
break;
case REGION2A:
fprintf(fp,"coordinates 2 A\n");
break;
case REGION1B:
fprintf(fp,"coordinates 1 B\n");
break;
case REGION2B:
fprintf(fp,"coordinates 2 B\n");
break;
}
/* do atoms */
for (list=model->cores ; list ; list=g_slist_next(list))
{
core = (struct core_pak *) list->data;
/* don't print checks */
if (core->region != region)
continue;
if (core->status & DELETED)
continue;
/* NB: x,y,z have been fractionalized */
ARR3SET(vec, core->x);
vecmat(model->latmat, vec);
/* NEW - compute molecule number */
if (core->bonds)
mol = g_slist_index(model->moles, core->mol);
else
mol = 0;
fprintf(fp,"%-7s core %14.9f %14.9f %14.9f %4d %14.9f\n",
core->atom_label, vec[0], vec[1], vec[2], mol, core->charge);
/* shell */
if (core->shell)
{
shel = (struct shel_pak *) core->shell;
ARR3SET(vec, shel->x);
vecmat(model->latmat, vec);
fprintf(fp,"%-7s shel %14.9f %14.9f %14.9f %4d %14.9f\n",
shel->shell_label, vec[0], vec[1], vec[2], mol, shel->charge);
}
}
/* done region */
fprintf(fp,"end\n");
}
/* done */
fclose(fp);
return(0);
}
/****************************/
/* marvin output file parse */
/****************************/
#define DEBUG_LOAD_MVNOUT 0
gint read_mvnout(gchar *filename, struct model_pak *model)
{
gboolean overwrite=FALSE;
gint i, j, region, offset, num_tokens;
gint count[REGION2B];
gchar *line, **buff, *text;
gdouble veca[2], vecb[2];
gdouble version;
gpointer scan;
GSList *core_list=NULL, *shell_list=NULL;
struct core_pak *core;
struct shel_pak *shell;
/* checks */
g_assert(model != NULL);
g_assert(filename != NULL);
scan = scan_new(filename);
if (!scan)
return(1);
/* setup */
model->id = MVNOUT;
model->periodic = 2;
strcpy(model->filename, filename);
g_free(model->basename);
model->basename = strdup_basename(filename);
for (i=REGION2B ; i-- ; )
count[i] = 0;
/* NB: disallow image creation in this direction */
model->pbc[2] = 0.0;
/* scan the file */
offset=0;
line = scan_get_line(scan);
while (!scan_complete(scan))
{
buff = tokenize(line, &num_tokens);
/* YAMH - get version to determine a parsing parameter */
if (g_ascii_strncasecmp(line, " Version", 27) == 0)
{
if (num_tokens > 1)
{
version = str_to_float(*(buff+1));
offset = (version < 1.99) ? 0 : 1;
#if DEBUG_LOAD_MVNOUT
printf("Marvin version %lf, offset = %d\n", version, offset);
#endif
}
}
/* get surface vectors */
if (g_ascii_strncasecmp(line,"Surface Vectors:",16) == 0)
{
/* read them in */
if (num_tokens > 3)
{
veca[0] = str_to_float(*(buff+2));
veca[1] = str_to_float(*(buff+3));
g_strfreev(buff);
buff = scan_get_tokens(scan, &num_tokens);
if (num_tokens > 1)
{
vecb[0] = str_to_float(*buff);
vecb[1] = str_to_float(*(buff+1));
/* NB: gdis stores vectors in reverse order */
VEC3SET(&(model->latmat[0]), veca[0], vecb[0], 0.0);
VEC3SET(&(model->latmat[3]), veca[1], vecb[1], 0.0);
VEC3SET(&(model->latmat[6]), 0.0, 0.0, 1.0);
model->construct_pbc = TRUE;
model->periodic = 2;
}
}
}
/* energy seaches */
if (g_ascii_strncasecmp(line, "Total Energy", 12) == 0)
{
g_strfreev(buff);
buff = g_strsplit(line, "=", 2);
text = format_value_and_units(*(buff+1), 5);
property_add_ranked(3, "Energy", text, model);
g_free(text);
}
if (g_ascii_strncasecmp(line, "Attachment Energy", 17) == 0)
{
if (num_tokens > 4)
{
if (model->gulp.eatt[0] == 0.0)
model->gulp.eatt[0] = str_to_float(*(buff+3));
else
model->gulp.eatt[1] = str_to_float(*(buff+3));
g_free(model->gulp.eatt_units);
model->gulp.eatt_units = g_strdup(*(buff+4));
}
}
if (g_ascii_strncasecmp(line, "Surface Energy", 14) == 0)
{
if (num_tokens > 5)
{
if (model->gulp.esurf[0] == 0.0)
model->gulp.esurf[0] = str_to_float(*(buff+3));
else
model->gulp.esurf[1] = str_to_float(*(buff+3));
g_free(model->gulp.esurf_units);
model->gulp.esurf_units = g_strdup(*(buff+5));
}
}
if (g_ascii_strncasecmp(line, "Gnorm", 5) == 0)
{
g_strfreev(buff);
buff = g_strsplit(line, "=", 2);
text = format_value_and_units(*(buff+1), 5);
property_add_ranked(6, "Gnorm", text, model);
g_free(text);
}
/* coordinate data search */
if (num_tokens > 5 && g_ascii_strncasecmp(line, "BLOCK", 5) == 0)
{
if (g_ascii_strncasecmp(*(buff+2),"region",5) == 0)
{
/* get region type */
region = REGION1A;
if (g_ascii_strncasecmp(*(buff+5),"1b",2) == 0)
region = REGION1B;
if (g_ascii_strncasecmp(*(buff+5),"2a",2) == 0)
region = REGION2A;
if (g_ascii_strncasecmp(*(buff+5),"2b",2) == 0)
region = REGION2B;
/* check if it's coordinate data (don't want gradients) */
line = scan_get_line(scan);
g_strfreev(buff);
buff = scan_get_tokens(scan, &num_tokens);
if (num_tokens > 3)
if (g_ascii_strncasecmp(*(buff+3),"coordinates",11) == 0)
{
model->region_empty[region] = FALSE;
/* only do this once as we don't want to overwrite from */
/* the start of the core/shell list at every new block */
if (count[region] && !overwrite)
{
#if DEBUG_LOAD_MVNOUT
printf("Repeat found - overwriting from now on.\n");
#endif
core_list = model->cores = g_slist_reverse(model->cores);
shell_list = model->shels = g_slist_reverse(model->shels);
overwrite = TRUE;
}
#if DEBUG_LOAD_MVNOUT
printf("Reading region: %d\n", region);
#endif
/* skip to data */
line = scan_get_line(scan);
line = scan_get_line(scan);
while (!scan_complete(scan))
{
if (g_ascii_strncasecmp(line,"-------",7) == 0)
break;
g_strfreev(buff);
buff = tokenize(line, &num_tokens);
if (num_tokens < 7)
break;
/* NEW - overwrite if 2nd occurance and existing, else create new */
/* NB: marvin 1.8 & 2.0 have this at different positions! */
if (g_ascii_strncasecmp(*(buff+offset+2),"cor",3) == 0)
{
if (overwrite && core_list)
{
core = core_list->data;
core_list = g_slist_next(core_list);
}
else
{
core = new_core(*(buff+offset+1), model);
model->cores = g_slist_prepend(model->cores, core);
core->charge = str_to_float(*(buff+offset+6));
core->lookup_charge = FALSE;
core->region = region;
}
core->x[0] = str_to_float(*(buff+offset+3));
core->x[1] = str_to_float(*(buff+offset+4));
core->x[2] = str_to_float(*(buff+offset+5));
}
if (g_ascii_strncasecmp(*(buff+offset+2),"she",3) == 0)
{
if (overwrite && shell_list)
{
shell = shell_list->data;
shell_list = g_slist_next(shell_list);
}
else
{
shell = new_shell(*(buff+offset+1), model);
model->shels = g_slist_prepend(model->shels, shell);
shell->charge = str_to_float(*(buff+offset+6));
shell->lookup_charge = FALSE;
shell->region = region;
}
shell->x[0] = str_to_float(*(buff+offset+3));
shell->x[1] = str_to_float(*(buff+offset+4));
shell->x[2] = str_to_float(*(buff+offset+5));
}
line = scan_get_line(scan);
}
count[region]++;
}
}
}
line = scan_get_line(scan);
g_strfreev(buff);
}
/* init lighting */
j=0;
for (i=REGION1A ; i<REGION2B ; i++)
{
if (!model->region_empty[i] && !j)
{
model->region[i] = TRUE;
j++;
}
else
model->region[i] = FALSE;
}
#if DEBUG_LOAD_MVNOUT
printf("Surface vectors: %lf x %lf\n",model->pbc[0],model->pbc[1]);
#endif
/* force non periodic in z (marvin is 2D) */
model->pbc[2] = 0;
model_prep(model);
scan_free(scan);
return(0);
}
/****************/
/* file reading */
/****************/
#define DEBUG_READ_MARVIN 0
gint read_marvin(gchar *filename, struct model_pak *model)
{
gint i, num_tokens, keyword, region;
gchar **buff, *tmp, line[LINELEN];
gdouble lat_mult = 1.0;
GSList *list;
struct core_pak *core;
struct shel_pak *shel;
FILE *fp;
/* checks */
g_return_val_if_fail(model != NULL, 1);
g_return_val_if_fail(filename != NULL, 1);
fp = fopen(filename, "rt");
if (!fp)
return(1);
/* defaults */
/*
model->periodic = 2;
*/
model->region_empty[REGION1A] = TRUE;
model->region_empty[REGION2A] = TRUE;
model->region_empty[REGION1B] = TRUE;
model->region_empty[REGION2B] = TRUE;
keyword = COMPARE;
region = REGION1A;
while(!fgetline(fp, line))
{
list = get_keywords(line);
if (list != NULL)
{
keyword = GPOINTER_TO_INT(list->data);
switch(keyword)
{
case MARVIN_COORD:
tmp = g_strstrip(g_strdup(get_token_pos(line, 1)));
if (g_ascii_strncasecmp("1 a", tmp, 3) == 0)
region = REGION1A;
if (g_ascii_strncasecmp("2 a", tmp, 3) == 0)
region = REGION2A;
if (g_ascii_strncasecmp("1 b", tmp, 3) == 0)
region = REGION1B;
if (g_ascii_strncasecmp("2 b", tmp, 3) == 0)
region = REGION2B;
model->region_empty[region] = FALSE;
g_free(tmp);
break;
}
}
else
{
buff = tokenize(line, &num_tokens);
/* have we got a valid element to read */
if (buff)
{
if (elem_symbol_test(*buff))
{
switch(keyword)
{
case MARVIN_COORD:
case MARVIN_BASIS:
if (num_tokens > 4)
{
if (g_ascii_strncasecmp("cor", *(buff+1), 3) == 0)
{
core = new_core(*buff, model);
model->cores = g_slist_prepend(model->cores, core);
core->region = region;
core->x[0] = str_to_float(*(buff+2));
core->x[1] = str_to_float(*(buff+3));
core->x[2] = str_to_float(*(buff+4));
if (num_tokens > 6)
{
core->charge = str_to_float(*(buff+6));
core->lookup_charge = FALSE;
}
}
if (g_ascii_strncasecmp("shel", *(buff+1), 4) == 0)
{
shel = new_shell(*buff, model);
model->shels = g_slist_prepend(model->shels, shel);
shel->region = region;
shel->x[0] = str_to_float(*(buff+2));
shel->x[1] = str_to_float(*(buff+3));
shel->x[2] = str_to_float(*(buff+4));
if (num_tokens > 6)
{
shel->charge = str_to_float(*(buff+6));
shel->lookup_charge = FALSE;
}
}
}
break;
}
}
else
{
/* Added by C. Fisher 2004 */
if (g_ascii_strncasecmp(line, "lattice", 7) == 0)
{
buff = tokenize(line, &num_tokens);
lat_mult = str_to_float(*(buff+1));
model->fractional = TRUE;
}
if (g_ascii_strncasecmp(line, "latvec", 6) == 0)
{
model->periodic = 2;
if (!fgetline(fp, line))
{
g_strfreev(buff);
buff = tokenize(line, &num_tokens);
model->latmat[0] = str_to_float(*(buff+0));
model->latmat[1] = str_to_float(*(buff+1));
model->latmat[2] = str_to_float(*(buff+2));
}
if (!fgetline(fp, line))
{
g_strfreev(buff);
buff = tokenize(line, &num_tokens);
model->latmat[3] = str_to_float(*(buff+0));
model->latmat[4] = str_to_float(*(buff+1));
model->latmat[5] = str_to_float(*(buff+2));
}
if (!fgetline(fp, line))
{
g_strfreev(buff);
buff = tokenize(line, &num_tokens);
model->latmat[6] = str_to_float(*(buff+0));
model->latmat[7] = str_to_float(*(buff+1));
model->latmat[8] = str_to_float(*(buff+2));
}
model->construct_pbc = TRUE;
}
/* Added by C. Fisher 2004 */
if (g_ascii_strncasecmp(line, "surfvec", 7) == 0)
{
if (!fgetline(fp, line))
{
g_strfreev(buff);
buff = tokenize(line, &num_tokens);
model->image_limit[1] = str_to_float(*(buff+0));
}
if (!fgetline(fp, line))
{
g_strfreev(buff);
buff = tokenize(line, &num_tokens);
model->image_limit[3] = str_to_float(*(buff+1));
}
}
}
}
g_strfreev(buff);
}
}
/* Added by C. Fisher 2004 */
for(i=0;i<9;i++)
model->latmat[i] *= lat_mult;
/* model setup */
strcpy(model->filename, filename);
g_free(model->basename);
model->basename = strdup_basename(filename);
model_prep(model);
/* Added by C. Fisher 2004 */
if( model->image_limit[1] != 1.0 || model->image_limit[3] != 1.0)
{
space_make_images(CREATE, model);
coords_init(CENT_COORDS, model);
}
fclose(fp);
return(0);
}
syntax highlighted by Code2HTML, v. 0.9.1