/* 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 */ /* * If you have any questions or comments about DL_POLY file support in GDIS, * send it to GDIS mailing list * http://lists.sourceforge.net/lists/listinfo/gdis-users * (you can also CC to Marcin Wojdyr, search the web for e-mail) * * Details about * The DL_POLY Molecular Simulation Package * can be found at http://www.cse.clrc.ac.uk/msi/software/DL_POLY/ * Following DL_POLY configuration files are (partially) supported: * * CONFIG/REVCON read-write * HISTORY formatted read * HISTORY unformatted (binary) not implemented * HISTORY file provides trajectory and can be very large. * * To see what features are not supported yet read comments and the code. * * When using shell model, there is no way to recognize which "atom" * represents core and which shell using only CONFIG/HISTORY file. * Convention is used: name of shell "atom" is created with * one of following postfixes: -shell _shell -shel _shel -shl _shl -sh _sh * eg. Zn-shl or Zn_shel */ #include #include #include #include "gdis.h" #include "coords.h" #include "model.h" #include "file.h" #include "parse.h" #include "scan.h" #include "matrix.h" #include "interface.h" /* main structures */ extern struct sysenv_pak sysenv; extern struct elem_pak elements[]; /* Replaces all blanks in str with '\0' and put pointers to tokens into tokens. * Returns number of tokens (n <= max_split). It doesn't allocates any memory. * tokens[] should have length >= max_split. */ int split_string(char *str, int max_split, char **tokens) { char *ptr; int counter=0; int was_space=1; for (ptr=str; *ptr; ++ptr) { if (isspace(*ptr)) { *ptr=0; was_space=1; } else if (was_space) { was_space=0; tokens[counter] = ptr; counter++; if (counter >= max_split) break; } } return counter; } int get_splitted_line(FILE *fp, int max_split, char **tokens) { gchar line[1024]; if (fgets(line, 1024, fp)) return split_string(line, max_split, tokens); else return 0; } void mark_animation_frame(struct model_pak *model, FILE *fp) { fpos_t offset; fgetpos(fp, &offset); model->frame_list = g_list_append(model->frame_list, g_memdup(&offset, sizeof(fpos_t))); model->num_frames++; } /**************** file reading ****************/ /* read boundary conditions */ static gint read_dlpoly_pbc(gint imcon, FILE *fp, struct model_pak *model) { /* imcon -- periodic boundary key: imcon meaning 0 no periodic boundaries 1 cubic boundary conditions 2 orthorhombic boundary conditions 3 parallelepiped boundary conditions 4 truncated octahedral boundary conditions 5 rhombic dodecahedral boundary conditions 6 x-y parallelogram boundary conditions with no periodicity in the z direction 7 hexagonal prism boundary conditions */ gint num_tokens, i, j; gchar *tokens[64]; gchar line[1024]; switch (imcon) { case 0: model->periodic = 0; return 0; case 1: case 2: model->periodic = 3; for (i=0; i < 3; ++i) { num_tokens = get_splitted_line(fp, 64, tokens); if (num_tokens < 3) { gui_text_show(ERROR, "not enought numbers in PBC desc.\n"); return 1; } for (j=0; j<3; ++j) model->latmat[3*i+j] = str_to_float(tokens[j]); } model->construct_pbc = TRUE; matrix_lattice_init(model); return 0; case 6: /*TODO*/ case 3: case 4: case 5: case 7: default: model->periodic = 0; /*eat 3 lines*/ for (i=0; i < 3; ++i) fgets(line, 1024, fp); gui_text_show(WARNING, "Sorry, not supported imcon (PBC key). " "Assuming no periodic boundaries.\n"); return 0; } } /* checks if "atom" is a shell in core-shell model */ static gint read_dlpoly_is_shell(gchar *name) { /* check for postfixes: -shell _shell -shel _shel -shl _shl -sh _sh * and return postfix length */ int len = strlen(name); if (len > 6 && (g_ascii_strcasecmp(name+len-6, "-shell") == 0 || g_ascii_strcasecmp(name+len-6, "_shell") == 0)) return 6; else if (len > 5 && (g_ascii_strcasecmp(name+len-5, "-shel") == 0 || g_ascii_strcasecmp(name+len-5, "_shel") == 0)) return 5; else if (len > 4 && (g_ascii_strcasecmp(name+len-4, "-shl") == 0 || g_ascii_strcasecmp(name+len-4, "_shl") == 0)) return 4; else if (len > 3 && (g_ascii_strcasecmp(name+len-3, "-sh") == 0 || g_ascii_strcasecmp(name+len-3, "_sh") == 0)) return 3; else return 0; } /* OH -> O, etc */ static void dlpoly_cut_elem_postfix(gchar *elem) { /* don't cut Cl, Ca, Cs, Cr, Co, Cu, Cd, Ce, Os, Hf */ gchar c = elem[0]; gchar x = 0; if (strlen(elem) <= 1) return; x = elem[1]; if (((c=='C' || c=='c') && (x!='l' && x!='L' && x!='a' && x!='A' && x!='s' && x!='S' && x!='r' && x!='R' && x!='o' && x!='O' && x!='u' && x!='U' && x!='d' && x!='D' && x!='e' && x!='E')) || ((c=='O' || c=='o') && (x!='s' && x!='S')) || ((c=='H' || c=='h') && (x!='f' && x!='F'))) elem[1] = '\0'; } static gint read_dlpoly_atoms(gint levcfg, FILE *fp, struct model_pak *model) { /* Possible values of levcfg: 1 - only coordinates in file 2 - coordinates and velocities 3 - coordinates, velocities and forces */ gchar elem[9]; gchar *tokens[64]; gint rec=0, num_tokens, sl=0; gdouble *x=0, *v=0; struct core_pak *core=NULL; struct shel_pak *shell=NULL; GSList *core_list, *shell_list; model->cores = g_slist_reverse(model->cores); model->shels = g_slist_reverse(model->shels); core_list = model->cores; shell_list = model->shels; while ((num_tokens = get_splitted_line(fp, 64, tokens))) { if (g_ascii_strncasecmp(tokens[0], "timestep", 8) == 0) { break; } switch (rec) { case 0: /* record i */ if (num_tokens < 1) { gui_text_show(ERROR, "Unexpected syntax in DL_POLY file" " (in record i).\n"); return 1; } sl = read_dlpoly_is_shell(tokens[0]); if (sl) { g_strlcpy(elem, tokens[0], 9); elem[strlen(elem) - sl] = '\0'; dlpoly_cut_elem_postfix(elem); if (shell_list) { shell = (struct shel_pak *) shell_list->data; shell_list = g_slist_next(shell_list); } else { shell = shell_new(elem, tokens[0], model); model->shels = g_slist_prepend(model->shels, shell); } x = shell->x; v = shell->v; } else { //core g_strlcpy(elem, tokens[0], 9); dlpoly_cut_elem_postfix(elem); if (core_list) { core = (struct core_pak *) core_list->data; core_list = g_slist_next(core_list); } else { core = core_new(elem, tokens[0], model); model->cores = g_slist_prepend(model->cores, core); } x = core->x; v = core->v; } break; case 1: /* record ii */ if (num_tokens < 3) { gui_text_show(ERROR, "Unexpected syntax in DL_POLY file" " (in record ii).\n"); return 1; } VEC4SET(x, str_to_float(tokens[0]), str_to_float(tokens[1]), str_to_float(tokens[2]), 1.0); break; case 2: /* record iii */ if (num_tokens < 3) { gui_text_show(ERROR, "Unexpected syntax in DL_POLY file" " (in record iii).\n"); return 1; } VEC3SET(v, str_to_float(tokens[0]), str_to_float(tokens[1]), str_to_float(tokens[2])); break; case 3: /* record iv */ if (num_tokens < 3) { gui_text_show(ERROR, "Unexpected syntax in DL_POLY file" " (in record iv).\n"); return 1; } /* no force info in GDIS */ break; default: /* never comes here */ printf("ERROR\n"); } rec = (rec+1) % (levcfg+2); /* rec = 0, 1, ..., levcfg+1, 0, 1, ... */ } model->cores = g_slist_reverse(model->cores); model->shels = g_slist_reverse(model->shels); printf("# cores: %i; # shells: %i\n", g_slist_length(model->cores), g_slist_length(model->shels)); return 0; } gint read_dlpoly(gchar *filename, struct model_pak *model) { gint num_tokens, history; gchar line[1024], *tokens[64]; FILE *fp; gint levcfg, /* tells if velocities and forces are in file*/ imcon; /* periodic boundary key */ int i; fp = fopen(filename, "rt"); if (!fp) return 1; /* first line is a title */ if (!fgets(line, 1024, fp)) return 2; /* cut trailing blanks */ for (i=strlen(line)-1; i>0 && isspace(line[i]); --i) line[i] = 0; model->title = g_strdup(line); /* second line */ num_tokens = get_splitted_line(fp, 64, tokens); if (num_tokens < 2) return 3; levcfg = (gint) str_to_float(tokens[0]); imcon = (gint) str_to_float(tokens[1]); /* 3rd line -- is this CONFIG or HISTORY? */ if (!fgets(line, 1024, fp)) return 4; history = (g_ascii_strncasecmp(line, "timestep", 8) == 0); if (history) { /* HISTORY */ /* read first frame */ num_tokens = split_string(line, 64, tokens); if (num_tokens >= 6) ; /* TODO read current nstep and time-step */ read_dlpoly_pbc(imcon, fp, model); read_dlpoly_atoms(levcfg, fp, model); /* mark all frames offsets */ model->num_frames = 0; fseek(fp, 0, SEEK_SET); /* TODO optimize it */ while (fgets(line, 1024, fp)) { if (g_ascii_strncasecmp(line, "timestep", 8) == 0) { fseek(fp, -strlen(line), SEEK_CUR); mark_animation_frame(model, fp); fseek(fp, strlen(line), SEEK_CUR); } } } else { /* CONFIG/REVCON */ if (num_tokens >= 4) ; /*TODO read current nstep and time-step */ read_dlpoly_pbc(imcon, fp, model); read_dlpoly_atoms(levcfg, fp, model); } /* model setup */ strcpy(model->filename, filename); g_free(model->basename); model->basename = strdup_basename(filename); model->fractional = FALSE; model_prep(model); return 0; } /******** animation -- overwrite frame ********/ gint read_dlpoly_frame(FILE *fp, struct model_pak *model) { gint num_tokens; gint imcon, levcfg; gchar *tokens[64]; puts("in read_dlpoly_frame()"); g_assert(fp != NULL); num_tokens = get_splitted_line(fp, 64, tokens); if (num_tokens < 6) return 11; g_assert(g_ascii_strncasecmp(tokens[0], "timestep", 8) == 0); levcfg = (gint) str_to_float(tokens[3]); imcon = (gint) str_to_float(tokens[4]); read_dlpoly_pbc(imcon, fp, model); read_dlpoly_atoms(levcfg, fp, model); return 0; } /**************** file writing ****************/ gint write_dlpoly(gchar *filename, struct model_pak *model) { GSList *list; struct core_pak *core; struct shel_pak *shell; FILE *fp; gint levcfg=0, imcon=0; gdouble x[3]; int i; /* 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); /* first line is a title */ fprintf(fp, "%s\n", model->title); if (model->periodic == 0) { imcon = 0; //no PBC } else if (model->periodic == 3) { if (model->latmat[0*3+1] == 0. && model->latmat[0*3+2] == 0. && model->latmat[1*3+0] == 0. && model->latmat[1*3+2] == 0. && model->latmat[2*3+0] == 0. && model->latmat[2*3+1] == 0.) { if (model->latmat[0*3+0] == model->latmat[1*3+1] && model->latmat[0*3+0] == model->latmat[2*3+2]) imcon = 1; //cubic else imcon = 2; //orthorhombic } /*TODO other cases */ } /* TODO when levcfg should be 1? */ fprintf(fp, "%10i%10i\n", levcfg, imcon); /* TODO write current nstep and time-step */ if (imcon > 0) { /*write PBC*/ for (i=0; i<3; ++i) fprintf(fp, "%20.10f%20.10f%20.10f\n", model->latmat[3*i], model->latmat[3*i+1], model->latmat[3*i+2]); } for (list=model->cores ; list ; list=g_slist_next(list)) { core = list->data; if (core->status & DELETED) continue; fprintf(fp, "%-8s\n", core->atom_label); ARR3SET(x, core->x); vecmat(model->latmat, x); fprintf(fp, "%20.8f%20.8f%20.8f\n", x[0], x[1], x[2]); if (levcfg > 0) fprintf(fp,"%20.8f%20.8f%20.8f\n",core->v[0],core->v[1],core->v[2]); } /*TODO order of atoms (shells and cores) should be preserved, * now if there are shells, they go at the end */ for (list=model->shels ; list ; list=g_slist_next(list)) { shell = list->data; if (shell->status & DELETED) continue; fprintf(fp, "%-8s\n", shell->shell_label); ARR3SET(x, shell->x); vecmat(model->latmat, x); fprintf(fp, "%20.8f%20.8f%20.8f\n", x[0], x[1], x[2]); if (levcfg > 0) fprintf(fp, "%20.8f%20.8f%20.8f\n", shell->v[0], shell->v[1], shell->v[2]); } fclose(fp); return 0; }