/*
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 <stdio.h>
#include <string.h>
#include "gdis.h"
#include "coords.h"
#include "file.h"
#include "parse.h"
#include "matrix.h"
#include "model.h"
#include "interface.h"
/* main structures */
extern struct sysenv_pak sysenv;
extern struct elem_pak elements[];
/* global variables */
gint last_frame;
/************************/
/* GAUSSIAN file saving */
/************************/
gint write_gauss(gchar *filename, struct model_pak *model)
{
gint i;
gdouble x[3];
GSList *list;
struct core_pak *core;
FILE *fp;
/* checks */
g_return_val_if_fail(model != NULL, 1);
g_return_val_if_fail(filename != NULL, 2);
/* open the file */
fp = fopen(filename,"wt");
if (!fp)
return(3);
g_free(model->basename);
model->basename = strdup_basename(filename);
fprintf(fp, "%%chk=%s.chk\n", model->basename);
fprintf(fp, "\n");
fprintf(fp, "# sp hf/sto-3g\n");
fprintf(fp, "\n");
gdis_blurb(fp);
if (model->title)
fprintf(fp, "%s\n", model->title);
fprintf(fp, "\n");
fprintf(fp, "0 1\n");
for (list=model->cores ; list ; list=g_slist_next(list))
{
core = (struct core_pak *) list->data;
if (core->status & DELETED)
continue;
/* everything is cartesian after latmat mult */
ARR3SET(x, core->x);
vecmat(model->latmat, x);
fprintf(fp, "%-7s %14.9f %14.9f %14.9f\n", elements[core->atom_code].symbol, x[0], x[1], x[2]);
}
/* write out lattice vectors */
if (model->periodic > 0)
for (i=0; i<model->periodic; i++)
fprintf(fp, "%-7s %14.9f %14.9f %14.9f\n", "Tv", model->latmat[0+i], model->latmat[3+i], model->latmat[6+i]);
fprintf(fp, "\n");
fclose(fp);
return(0);
}
gchar **strip_empty_strings(gchar **str_array)
{
gchar **new_str_array = NULL;
gint i=0, j=0;
while (str_array[i] != NULL)
{
if (strlen(str_array[i]) != 0)
{
new_str_array = g_realloc(new_str_array, (j + 2)*sizeof(char *)); /* add one for NULL pointer */
new_str_array[j] = str_array[i];
j++;
}
else
g_free(str_array[i]);
i++;
}
new_str_array[j] = NULL;
g_free(str_array);
return(new_str_array);
}
gchar **read_gauss_keyword(gchar *string)
{
gint i;
/*
Options to keywords may be specified in any of the following forms:
ÊÊÊÊÊÊÊÊÊ keyword = option
ÊÊÊÊÊÊÊÊÊÊÊÊkeyword(option)
ÊÊÊÊÊÊÊÊÊÊÊkeyword=(option1, option2, ...)
ÊÊÊÊÊÊÊÊÊÊÊÊkeyword(option1, option2, ...)
*/
/* Spaces, tabs, commas, or forward slashes can be used in any combination to separate items within a line. Multiple spaces are treated as a single delimiter. */
for (i=0; i<strlen(string); i++)
switch (string[i])
{
case '=':
case '(':
case ')':
case '\t':
case ',': string[i] = ' ';
}
return(strip_empty_strings(g_strsplit(string, " ", 0)));
}
/*******************************************/
/* read single GAUSSIAN output configuration */
/*******************************************/
#define DEBUG_READ_GAUSSIAN_OUT 0
gint read_gauss_out_block(FILE *fp, struct model_pak *model)
{
gint i;
gint num_tokens;
gchar **buff, line[LINELEN], *text;
GSList *clist;
GString *energy_string, *grad_string;
struct core_pak *core;
clist = model->cores;
model->periodic = 0;
model->fractional = FALSE;
model->construct_pbc = FALSE;
/* 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, " GradGradGradGrad", 16) == 0)
break;*/
/* read coordinates */
if (g_strrstr(line, "Center Atomic") != NULL)
{
for (i=0; i<3; i++)
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 == 6)
{
/* atomic number greater than 0 => real atom */
if (str_to_float(*(buff+1)) > 0)
{
if (clist)
{
core = (struct core_pak *) clist->data;
clist = g_slist_next(clist);
}
else
{
core = new_core(elements[(gint) str_to_float(*(buff+1))].symbol, 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_GAUSSIAN_OUT
printf("new coords %f %f %f\n", core->x[0], core->x[1], core->x[2]);
#endif
}
/* atomic number = -2 => lattice vector */
if (str_to_float(*(buff+1)) == -2)
{
model->latmat[0+model->periodic] = str_to_float(*(buff+3));
model->latmat[3+model->periodic] = str_to_float(*(buff+4));
model->latmat[6+model->periodic] = str_to_float(*(buff+5));
model->periodic++;
model->construct_pbc = TRUE;
#if DEBUG_READ_GAUSSIAN_OUT
printf(">> latmat: ");
P3MAT(" ", model->latmat);
#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;
}
/* Energy */
if (g_strrstr(line, "SCF Done:") != NULL)
{
buff = g_strsplit(line, "=", 2);
text = format_value_and_units(g_ascii_strdown(*(buff+1), -1), 7);
property_add_ranked(3, "Energy", text, model);
g_free(text);
g_strfreev(buff);
}
if (g_strrstr(line, "EUMP2") != NULL)
{
buff = g_strsplit(line, "=", 3);
energy_string = g_string_new("");
g_string_append_printf(energy_string, "%.7f %s", str_to_float(*(buff+2)), "a.u.");
property_add_ranked(3, "Energy", energy_string->str, model);
g_string_free(energy_string, TRUE);
g_strfreev(buff);
}
/* TODO higher order MP */
/* gradients */
if (g_strrstr(line, "Internal Forces:") != NULL)
{
buff = tokenize(line, &num_tokens);
grad_string = g_string_new("");
g_string_append_printf(grad_string, "%.7f %s", str_to_float(*(buff+5)), "a.u./Bohr");
property_add_ranked(4, "RMS Gradient", grad_string->str, model);
g_string_free(grad_string, TRUE);
g_strfreev(buff);
}
if (g_strrstr(line, "Predicted change in Energy") != NULL)
{
if (fgetline(fp, line))
{
gui_text_show(ERROR, "unexpected end of file checking for succesful optimization\n");
return(2);
}
if (g_ascii_strcasecmp(property_lookup("Calculation", model), "Forces") == 0)
last_frame = TRUE;
if (g_strrstr(line, "Optimization completed") != NULL)
last_frame = TRUE;
/* fix me - better clue to stopping */
break;
}
if (g_strrstr(line, "Job cpu time") != NULL)
{
last_frame = TRUE;
/* fix me - better clue to stopping */
break;
}
}
return(0);
}
/*******************************/
/* GAUSSIAN output frame reading */
/*******************************/
gint read_gauss_out_frame(FILE *fp, struct model_pak *model)
{
return(read_gauss_out_block(fp, model));
}
/********************************/
/* Read in a GAUSSIAN output file */
/********************************/
gint read_gauss_out(gchar *filename, struct model_pak *model)
{
gint frame, i, num_tokens;
gint have_scan=FALSE;
gchar **buff, **buff1, **tokens, line[LINELEN];
FILE *fp;
fp = fopen(filename, "rt");
if (!fp)
return(1);
last_frame = FALSE;
frame=0;
/* defaults */
property_add_ranked(2, "Calculation", "Single Point", model);
property_add_ranked(6, "Method", "HF", model);
property_add_ranked(7, "Basis", "STO-3G", model);
while (!fgetline(fp, line))
{
/* find gaussian route section and store selected parameters in hash table */
if (g_ascii_strncasecmp(line, " #", 2) == 0)
{
/* Spaces, tabs, commas, or forward slashes can be used in any combination to separate items within a line. Multiple spaces are treated as a single delimiter. */
buff = tokenize((gchar *)(line) + 2, &num_tokens);
for (i=0; i<num_tokens; i++)
{
/* handle iops */
if (g_strrstr(buff[i], "iop") != NULL)
continue;
/* method and basis handled specially */
/* hack to find method and basis but not general as one or other can be missing */
if (g_strrstr(buff[i], "/") != NULL)
{
buff1 = g_strsplit(buff[i], "/", 2);
property_add_ranked(6, "Method", buff1[0], model);
property_add_ranked(7, "Basis", buff1[1], model);
g_strfreev(buff1);
continue;
}
tokens = read_gauss_keyword(buff[i]);
if (g_ascii_strncasecmp(tokens[0], "Opt", 3) == 0)
{
property_add_ranked(2, "Calculation", "Optimisation", model);
}
else if (g_ascii_strncasecmp(tokens[0], "Force", 5) == 0)
{
property_add_ranked(2, "Calculation", "Forces", model);
continue;
}
else if (g_ascii_strncasecmp(tokens[0], "SCRF", 4) == 0)
{
property_add_ranked(6, "Solvation", "", model);
}
g_strfreev(tokens);
}
g_strfreev(buff);
}
/* ModRedundant Section */
if (g_ascii_strncasecmp(line, " The following ModRedundant input section has been read:", 56) == 0)
{
while (TRUE)
{
if (fgetline(fp, line))
{
gui_text_show(ERROR, "unexpected end of file reading ModRedundant input section\n");
return(2);
}
if (g_ascii_strncasecmp(line, " Iteration", 10) == 0)
break;
buff = tokenize(line, &num_tokens);
for (i=0; i<num_tokens; i++)
{
if (g_ascii_strncasecmp(buff[i], "S", 1) == 0)
have_scan = TRUE;
}
g_strfreev(buff);
}
}
/* read coordinates */
if ((g_strrstr(line, "Input orientation:") != NULL) || (g_strrstr(line, "Z-Matrix orientation:") != NULL) || (g_strrstr(line, "Standard orientation:") != NULL))
{
/* go through all frames to count them */
add_frame_offset(fp, model);
read_gauss_out_block(fp, model);
frame++;
if (!have_scan && last_frame)
break;
}
}
/* read extra data from repeat of last coordinates here */
/* 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