/* 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 #include #include #include #include "gdis.h" #include "coords.h" #include "error.h" #include "file.h" #include "parse.h" #include "matrix.h" #include "zmatrix.h" #include "zmatrix_pak.h" #include "model.h" #include "interface.h" gint set_true_false(); gint set_energy_units(); gint set_temperature_units(); gint set_length_units(); gint set_force_units(); gint set_pressure_units(); gint set_time_units(); gint set_mominert_units(); gint print_file_energy_units(); gint print_file_mass_units(); gint print_file_length_units(); gint print_file_time_units(); gint print_file_temperature_units(); gint print_file_pressure_units(); gint print_file_force_units(); gint print_file_mominert_units(); enum {SIESTA_DEFAULT, SIESTA_SPECIES, SIESTA_LATTICE_PARAMETERS, SIESTA_LATTICE_VECTORS, SIESTA_COORDS, SIESTA_ZMATRIX}; /* main structures */ extern struct sysenv_pak sysenv; extern struct elem_pak elements[]; /**********************************/ /* obtain a core's species number */ /**********************************/ gint fdf_species_index(gchar *atom_label, GSList *species_list) { gint i; struct species_pak *species; GSList *list; /* find corresponding number of this element in the species block */ i=1; for (list=species_list ; list ; list=g_slist_next(list)) { species = list->data; if (g_ascii_strcasecmp(atom_label, species->label) == 0) return(i); i++; } return(0); } /*****************************************************/ /* generate the list of unique species, SIESTA style */ /*****************************************************/ GSList *fdf_species_build(struct model_pak *model) { gint code; struct species_pak *species_data; GSList *list, *normal_list, *ghost_list, *species_list; /* init the normal, ghost, and overall species lists */ normal_list = find_unique(LABEL_NORMAL, model); ghost_list = find_unique(LABEL_GHOST, model); species_list = NULL; for (list=normal_list ; list ; list=g_slist_next(list)) { code = elem_symbol_test(list->data); /* update the species identifier list */ species_data = g_malloc(sizeof(struct species_pak)); species_data->label = list->data; species_data->number = code; species_list = g_slist_prepend(species_list, species_data); } for (list=ghost_list ; list ; list=g_slist_next(list)) { code = elem_symbol_test(list->data); /* update the species identifier list */ species_data = g_malloc(sizeof(struct species_pak)); species_data->label = list->data; species_data->number = -code; species_list = g_slist_prepend(species_list, species_data); } species_list = g_slist_reverse(species_list); /* parse the zmatrix entries and make the type match the (possible new) species order */ if (model->zmatrix) { struct zmat_pak *zmatrix = model->zmatrix; /* scan zval entries - make sure type value reflects element order */ for (list=zmatrix->zlines ; list ; list=g_slist_next(list)) { struct zval_pak *zval = list->data; zval->type = fdf_species_index(zval->elem, species_list); } } return(species_list); } /****************/ /* file writing */ /****************/ gint write_fdf(gchar *filename, struct model_pak *model) { gint i; gdouble x[3], depth; GSList *list, *clist, *species_list; struct core_pak *core; struct species_pak *species_data; 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); /* print header */ fprintf(fp, "# "); gdis_blurb(fp); fprintf(fp, "#\n\n"); fprintf(fp, "SystemLabel %s\n\n", model->basename); fprintf(fp, "NumberOfAtoms %d\n\n", g_slist_length(model->cores)); /* init the normal, ghost, and overall species lists */ species_list = fdf_species_build(model); fprintf(fp, "NumberOfSpecies %d\n", g_slist_length(species_list)); /* write the species (unique atom types) block */ fprintf(fp, "%%block ChemicalSpeciesLabel\n"); i=1; for (list=species_list ; list ; list=g_slist_next(list)) { species_data = list->data; fprintf(fp, " %3d %3d %s\n", i, species_data->number, species_data->label); i++; } fprintf(fp, "%%endblock ChemicalSpeciesLabel\n\n"); /* write the lattice data */ if (model->periodic) { fprintf(fp, "LatticeConstant 1.0 Ang\n"); if (model->construct_pbc) { /* NB: siesta matrices are transposed wrt gdis */ fprintf(fp, "%%block LatticeVectors\n"); fprintf(fp,"%15.10f %15.10f %15.10f\n", model->latmat[0], model->latmat[3], model->latmat[6]); fprintf(fp,"%15.10f %15.10f %15.10f\n", model->latmat[1], model->latmat[4], model->latmat[7]); fprintf(fp,"%15.10f %15.10f %15.10f\n", model->latmat[2], model->latmat[5], model->latmat[8]); fprintf(fp, "%%endblock LatticeVectors\n\n"); } else { fprintf(fp, "%%block LatticeParameters\n"); /* if it's a surface - print Dhkl (mult by region sizes) */ if (model->periodic == 2) { /* get depth info */ depth = (model->surface.region[0]+model->surface.region[1]) *model->surface.depth; /* no depth info - make it large enough to fit everything */ if (depth < POSITION_TOLERANCE) depth = 2.0*model->rmax; fprintf(fp, " %f %f %f", model->pbc[0], model->pbc[1], depth); } else fprintf(fp, " %f %f %f", model->pbc[0], model->pbc[1], model->pbc[2]); fprintf(fp, " %f %f %f\n", R2D*model->pbc[3], R2D*model->pbc[4], R2D*model->pbc[5]); fprintf(fp, "%%endblock LatticeParameters\n\n"); } } /* write the zmatrix data */ /* if (model->zmatrix) */ if (zmat_entries_get(model->zmatrix)) { fprintf(fp, "ZM.UnitsLength %s\n", zmat_distance_units_get(model->zmatrix)); fprintf(fp, "ZM.UnitsAngle %s\n\n", zmat_angle_units_get(model->zmatrix)); fprintf(fp, "%%block Zmatrix\n"); zmat_coord_print(fp, species_list, model); zmat_mol_print(fp, model->zmatrix); zmat_var_print(fp, model->zmatrix); zmat_const_print(fp, model->zmatrix); fprintf(fp, "%%endblock Zmatrix\n"); } else { /* write the atoms - for surfaces (and isolated molecules), coords must be cartesian */ if (model->fractional) fprintf(fp, "AtomicCoordinatesFormat Fractional\n"); else fprintf(fp, "AtomicCoordinatesFormat NotScaledCartesianAng\n"); fprintf(fp, "%%block AtomicCoordinatesAndAtomicSpecies\n"); for (clist=model->cores ; clist ; clist=g_slist_next(clist)) { core = clist->data; if (core->status & DELETED) continue; ARR3SET(x, core->x); /* NB: want fractional if 3D periodic, otherwise cartesian */ if (!model->fractional) vecmat(model->latmat, x); /* find corresponding number of this element in the species block */ i=1; for (list=species_list ; list ; list=g_slist_next(list)) { species_data = list->data; if (g_ascii_strcasecmp(core->atom_label, species_data->label) == 0) { if (core->ghost) { if (species_data->number == -core->atom_code) break; } else { if (species_data->number == core->atom_code) break; } } i++; } /* found? */ if (list) fprintf(fp," %14.9f %14.9f %14.9f %d\n", x[0], x[1], x[2], i); else printf("write_fdf() error: bad species type.\n"); } fprintf(fp, "%%endblock AtomicCoordinatesAndAtomicSpecies\n\n"); } free_slist(species_list); /* * Siesta option writer */ // PAO.BasisType switch (model->siesta.basis_type) { case SPLIT_PAOBT: { fprintf(fp, "PAO.BasisType split\n"); break; } case SPLITGAUSS_PAOBT: { fprintf(fp, "PAO.BasisType splitgauss\n"); break; } case NODES_PAOBT: { fprintf(fp, "PAO.BasisType nodes\n"); break; } case NONODES_PAOBT: { fprintf(fp, "PAO.BasisType nonodes\n"); break; } default: { //split! fprintf(fp, "PAO.BasisType split\n"); break; } } switch (model->siesta.basis_set) { case SZ_ZETA : { fprintf(fp,"PAO.BasisSize SZ\n"); break; } case SZP_ZETA : { fprintf(fp,"PAO.BasisSize SZP\n"); break; } case DZ_ZETA : { fprintf(fp,"PAO.BasisSize DZ\n"); break; } case DZP_ZETA : { fprintf(fp,"PAO.BasisSize DZP\n"); break; } case CUSTOM_ZETA : { //TODO - CUSTOM Zeta Levels? - should set a tag somewhere fprintf(fp,"PAO.BasisSize DZP\n"); break; } } //siesta_solution_method_type switch (model->siesta.solution_method_type) { case DIAGON_SOLMETH: fprintf(fp, "SolutionMethod diagon\n"); break; case ORDERN_SOLMETH: fprintf(fp, "SolutionMethod OrderN\n"); break; default: //should never get here.... error_table_entry("Inconsistancy model->siesta.solution_method_type.\n"); fprintf(fp, "SolutionMethod diagon\n"); break; } // if (model->siesta.split_zeta_norm != 0.15) fprintf(fp, "PAO.SplitNorm %f\n", model->siesta.split_zeta_norm); fprintf(fp, "PAO.EnergyShift %f ", model->siesta.energy_shift); if (model->siesta.energy_shift_units) print_file_energy_units(fp, model->siesta.energy_shift_units); else fprintf(fp, "\n"); //Default is FALSE. if (model->siesta.harris_functional == TRUE) fprintf(fp, "Harris_functional true\n"); else fprintf(fp, "Harris_functional false\n"); //Default == LDA switch (model->siesta.xc_functional_type) { case GGA_XCFUNC: fprintf(fp, "XC.functional GGA\n"); break; case LDA_XCFUNC: fprintf(fp, "XC.functional LDA\n"); break; case LSD_XCFUNC: fprintf(fp, "XC.functional LSD\n"); break; default: //should never get here error_table_entry("Inconsistancy model->siesta.xc_functional.\n"); break; } switch (model->siesta.xc_author_type) { case CA_XCAUTH: fprintf(fp, "XC.Authors CA\n"); break; case PZ_XCAUTH: fprintf(fp, "XC.Authors PZ\n"); break; case PW_XCAUTH: fprintf(fp, "XC.Authors PW92\n"); break; case PBE_XCAUTH: fprintf(fp, "XC.Authors PBE\n"); break; case RPBE_XCAUTH: fprintf(fp, "XC.Authors RPBE\n"); break; case LYP_XCAUTH: fprintf(fp, "XC.Authors LYP\n"); break; default: //should never get here error_table_entry("Inconsistancy model->siesta.xc_authors.\n"); break; } if (model->siesta.spin_polarised == TRUE) fprintf (fp, "SpinPolarized true\n"); else fprintf (fp, "SpinPolarized false\n"); //MeshCutoff fprintf (fp, "MeshCutoff %f", model->siesta.mesh_cutoff); if (print_file_energy_units(fp, model->siesta.mesh_cutoff_units) != 0) { error_table_entry("Warning - Units missing MeshCutoff\n"); } //kgrid_cutoff fprintf(fp, "kgrid_cutoff %f", model->siesta.kgrid_cutoff); if (print_file_length_units(fp, model->siesta.kgrid_cutoff_units) !=0) { error_table_entry("Warning - Units missing kgrid_cutoff\n"); } fprintf(fp, "ElectronicTemperature %f", model->siesta.electronic_temperature); if (model->siesta.electronic_temperature_units != -1) { if (print_file_temperature_units(fp, model->siesta.electronic_temperature_units) !=0) { error_table_entry("Warning - Units missing electronic_temperature_units\n"); } } fprintf(fp, "MaxSCFIterations %d\n", (int) model->siesta.no_of_cycles); fprintf(fp, "DM.NumberPulay %d\n", (int) model->siesta.no_of_pulay_matrices); fprintf(fp, "DM.MixingWeight %f\n", model->siesta.mixing_weight); if (model->siesta.pulay_mixing) { //what does this map to? } switch (model->siesta.run_type) { case SINGLE_POINT: fprintf(fp, "MD.TypeOfRun CG\n"); break; case OPTIMISATION: fprintf(fp, "MD.TypeOfRun Verlet\n"); break; case MOLECULAR_DYNAMICS: fprintf(fp, "MD.TypeOfRun Verlet\n"); break; case PHONON_CALCULATION: switch (model->siesta.md_type_of_run) { case FC_MDRUN: fprintf(fp, "MD.TypeOfRun FC\n"); break; case PHONON_MDRUN: fprintf(fp, "MD.TypeOfRun Phonon\n"); break; default: break; } break; //MORE OPTIONS HERE > see manual default : fprintf(fp," \n"); error_table_entry("Warning - MD.TypeOfRun *NOT* written out\n"); break; } if (model->siesta.md_variable_cell) fprintf(fp, "MD.VariableCell true\n"); else fprintf(fp, "MD.VariableCell false\n"); if (model->siesta.md_num_cg_steps != -1.0) { fprintf(fp, "MD.NumCGsteps %d\n", (gint) model->siesta.md_num_cg_steps); } fprintf(fp, "MD.MaxCGDispl %f", model->siesta.md_max_cg_displacement); if (print_file_length_units(fp, model->siesta.md_max_cg_displacement_units) != 0) { error_table_entry("Warning - Units missing MD.MaxCGDispl\n"); } fprintf(fp, "MD.PreconditionVariableCell %f", model->siesta.md_precondition_variable_cell); if (print_file_length_units(fp, model->siesta.md_precondition_variable_cell_units ) != 0) { error_table_entry("Warning - Units missing MD.PreconditionVariableCell\n"); } fprintf(fp, "MD.MaxStressTol %f", model->siesta.md_max_stress_tol); if (print_file_pressure_units(fp, model->siesta.md_max_stress_tol_units) != 0) { error_table_entry("Warning - Units missing MD.MaxStressTol\n"); } fprintf(fp, "\n%%block MD.TargetStress\n"); fprintf(fp, "%f %f %f %f %f %f\n", model->siesta.md_target_stress_xx, model->siesta.md_target_stress_yy, model->siesta.md_target_stress_zz, model->siesta.md_target_stress_xy, model->siesta.md_target_stress_xz, model->siesta.md_target_stress_yz); fprintf(fp, "%%endblock MD.TargetStress\n\n"); fprintf(fp, "MD.MaxForceTol %f", model->siesta.md_max_force_tol); if (print_file_force_units(fp, model->siesta.md_max_force_tol_units) != 0) { error_table_entry("Warning - Units missing MD.MaxForceTol\n"); } //run_type is not CG_MDRUN fprintf(fp, "MD.InitialTimeStep %d\n", (gint) model->siesta.md_inital_time_step); fprintf(fp, "MD.FinalTimeStep %d\n", (gint) model->siesta.md_final_time_step); fprintf(fp, "MD.LengthTimeStep %d", (gint) model->siesta.timestep); if (print_file_time_units(fp, model->siesta.timestep_units) != 0) { error_table_entry("Warning - Units missing MD.LengthTimeStep\n"); } fprintf(fp, "MD.InitialTemperature %f", model->siesta.md_inital_temperature); if (print_file_temperature_units(fp, model->siesta.md_inital_temperature_units) != 0) { error_table_entry("Warning - Units missing MD.InitialTemperature\n"); } if (model->siesta.md_quench) fprintf(fp, "MD.Quench true\n"); else fprintf(fp, "MD.Quench false\n"); fprintf(fp, "MD.TargetTemperature %f", model->siesta.md_inital_temperature); if (print_file_temperature_units(fp, model->siesta.md_inital_temperature_units) !=0) { error_table_entry("Warning - Units missing MD.TargetTemperature\n"); } fprintf(fp, "MD.NoseMass %f", model->siesta.md_nose_mass); if (print_file_mominert_units(fp, model->siesta.md_nose_mass_units) != 0) { error_table_entry("Warning - Units missing MD.NoseMass\n"); } fprintf(fp, "MD.ParrinelloRahmanMass %f", model->siesta.md_parrinello_rahman_mass); if (print_file_mominert_units(fp, model->siesta.md_parrinello_rahman_mass_units) != 0) { error_table_entry("Warning - Units missing MD.ParrinelloRahmanMass\n"); } switch (model->siesta.md_anneal_option) { case TEMPERATURE_ANNEAL: fprintf(fp, "MD.AnnealOption Temperature\n"); break; case PRESSURE_ANNEAL: fprintf(fp, "MD.AnnealOption Pressure\n"); break; case TEMP_AND_PRESS_ANNEAL: fprintf(fp, "MD.AnnealOption TemperatureandPressure\n"); break; default: fprintf(fp, "MD.AnnealOption TemperatureandPressure\n"); break; } fprintf(fp, "MD.TauRelax %f", model->siesta.md_tau_relax); if (print_file_time_units(fp, model->siesta.md_tau_relax_units) !=0) { error_table_entry("Warning - Units missing MD.TauRelax\n"); } fprintf(fp, "MD.BulkModulus %f", model->siesta.md_bulk_modulus); if (print_file_pressure_units(fp, model->siesta.md_bulk_modulus_units) !=0) { error_table_entry("Warning - Units missing MD.BulkModulus\n"); } fprintf(fp, "MD.TargetPressure %f", model->siesta.md_target_pressure); if (print_file_pressure_units(fp, model->siesta.md_target_pressure_units) !=0) { error_table_entry("Warning - Units missing MD.BulkModulus\n"); } fprintf(fp, "MD.FCDispl %f", model->siesta.md_fc_displ); if (print_file_length_units(fp, model->siesta.md_fc_displ_units) !=0) { error_table_entry("Warning - Units missing MD.FCDispl\n"); } fprintf(fp, "MD.FCfirst %d\n", model->siesta.md_fc_first ); fprintf(fp, "MD.FClast %d\n\n", model->siesta.md_fc_last); if (model->siesta.use_saved_data) fprintf(fp, "UseSaveData\t\ttrue\n"); else fprintf(fp, "UseSaveData\t\tfalse\n"); if (model->siesta.file_output_write_coor_init) //WriteCoorInital fprintf(fp, "WriteCoorInital true\n"); else fprintf(fp, "WriteCoorInital false\n"); if (model->siesta.file_output_write_coor_step) //WriteCoorStep fprintf(fp, "WriteCoorStep true\n"); else fprintf(fp, "WriteCoorStep false\n"); if (model->siesta.file_output_write_forces) //WriteForces fprintf(fp, "WriteForces\t\ttrue\n"); else fprintf(fp, "WriteForces\t\tfalse\n"); if (model->siesta.file_output_write_kpoints) //WriteKpoints fprintf(fp, "WriteKpoints\t\ttrue\n"); else fprintf(fp, "WriteKpoints\t\tfalse\n"); if (model->siesta.file_output_write_eigenvalues) //WriteEigenvalues fprintf(fp, "WriteEigenvalues\t\ttrue\n"); else fprintf(fp, "WriteEigenvalues\t\tfalse\n"); if (model->siesta.file_output_write_kbands) //WriteKbands fprintf(fp, "WriteKbands\t\ttrue\n"); else fprintf(fp, "WriteKbands\t\tfalse\n"); if (model->siesta.file_output_write_bands) //WriteBands fprintf(fp, "WriteBands\t\ttrue\n"); else fprintf(fp, "WriteBands\t\tfalse\n"); if (model->siesta.file_output_write_wavefunctions) //WriteWaveFunctions fprintf(fp, "WriteWaveFunctions\t\ttrue\n"); else fprintf(fp, "WriteWaveFunctions\t\tfalse\n"); if (model->siesta.file_output_write_mullikenpop != -1) //WriteMullikenPop fprintf(fp, "WriteMullikenPop\t\t%d\n", (int) model->siesta.file_output_write_mullikenpop); if (model->siesta.file_output_write_dm) //WriteDM fprintf(fp, "WriteDM\t\ttrue\n"); else fprintf(fp, "WriteDM\t\tfalse\n"); if (model->siesta.file_output_write_coor_xmol) //WriteCoorXmol fprintf(fp, "WriteCoorXmol\t\ttrue\n"); else fprintf(fp, "WriteCoorXmol\t\tfalse\n"); if (model->siesta.file_output_write_coor_cerius) //WriteCoorCerius fprintf(fp, "WriteCoorCerius\t\ttrue\n"); else fprintf(fp, "WriteCoorCerius\t\tfalse\n"); if (model->siesta.file_output_write_md_xmol) //WriteMDXmol fprintf(fp, "WriteMDXmol\t\ttrue\n"); else fprintf(fp, "WriteMDXmol\t\tfalse\n"); if (model->siesta.file_output_write_md_history) //WriteMDhistory fprintf(fp, "WriteMDhistory\t\ttrue\n"); else fprintf(fp, "WriteMDhistory\t\tfalse\n"); /* g_slist_free(normal_list); g_slist_free(ghost_list); */ fclose(fp); return(0); } /****************************************************/ /* attempt to read header at a particular precision */ /****************************************************/ gint siesta_header_read(FILE *fp, gint precision, gdouble *cell, gint *grid) { gint i, swap=FALSE; int grid_int[4]; /*long grid_long[4];*/ /*float cell_float[9];*/ double cell_double[9]; /* standard return values */ for (i=9 ; i-- ; ) cell[i] = 0.0; for (i=4 ; i-- ; ) grid[i] = 0; /* precision dependent read */ switch (precision) { case 1: READ_RECORD; fread(cell_double, sizeof(double), 9, fp); READ_RECORD; READ_RECORD; fread(&grid_int, sizeof(int), 4, fp); READ_RECORD; /* byte ordering check */ for (i=9 ; i-- ; ) { if (fabs(cell_double[i]) > 1000000.0) { swap = TRUE; break; } } if (swap) { for (i=9 ; i-- ; ) swap_bytes(&cell_double[i], sizeof(double)); for (i=4 ; i-- ; ) swap_bytes(&grid_int[i], sizeof(int)); } /* pass back */ for (i=9 ; i-- ; ) cell[i] = cell_double[i]; for (i=4 ; i-- ; ) grid[i] = grid_int[i]; break; default: printf("FIXME - unsuported precision.\n"); } return(swap); } /******************************************/ /* check for additional siesta data files */ /******************************************/ #define DEBUG_SIESTA_DENSITY 0 void seista_density_files(const gchar *source, struct model_pak *model) { gint grid[4]; gint i, j, k, n, size, file_size, type, swap; gchar *text, *name; float *f; gdouble cell[9], icell[9], *ptr; FILE *fp; g_assert(source != NULL); g_assert(model != NULL); #if DEBUG_SIESTA_DENSITY printf("source: %s\n", source); #endif /* search for siesta output types */ /* in same directory & with same basename - different extension */ for (type=0 ; type<2 ; type++) { /* remove extension */ name = text = NULL; n = strlen(source); for (i=n ; i-- ; ) { if (*(source+i) == '.') { text = g_strndup(source, i); break; } } /* should never get a name with no extension */ if (!text) return; /* build current type name */ switch (type) { case 0: name = g_strdup_printf("%s.RHO", text); break; case 1: name = g_strdup_printf("%s.VH", text); break; } g_free(text); if (!name) continue; fp = fopen(name, "rt"); if (!fp) { #if DEBUG_SIESTA_DENSITY printf("Could not find: %s\n", name); #endif continue; } file_size = file_byte_size(name); g_free(name); swap = siesta_header_read(fp, 1, cell, grid); memcpy(icell, cell, 9*sizeof(gdouble)); matrix_invert(icell); #if DEBUG_SIESTA_DENSITY printf("swap bytes: %d\n", swap); P3MAT(" cell: ", cell); P3MAT("icell: ", icell); printf("grid: %d x %d x %d : %d\n", grid[0], grid[1], grid[2], grid[3]); #endif /* NB: FORTRAN has a 4 byte separator either side of each record (single write() call) */ /* real*8 cell(3,3) real*4 rho(nmesh(1)*nmesh(2)*nmesh(3),nspin) integer mesh(3) integer nspin write(iu) cell write(iu) mesh, nspin ind = 0 do is = 1,nspin do iz = 1,mesh(3) do iy = 1,mesh(2) write(iu) rho(ind+ix,is),ix = 1,mesh(1) ind = ind + mesh(1) enddo enddo enddo */ /* data array size */ n = grid[1] * grid[2] * grid[3]; size = 8 + 4*grid[0]; size *= n; /* cell + mesh size */ size += 104; #if DEBUG_SIESTA_DENSITY printf("Data size: %d bytes\n", size); printf("File size: %d bytes\n", file_size); #endif if (size != file_size) { /* TODO - retry at different precision */ gui_text_show(ERROR, "Failed to parse binary file.\n"); continue; } memcpy(model->siesta.cell, cell, 9*sizeof(gdouble)); memcpy(model->siesta.icell, icell, 9*sizeof(gdouble)); memcpy(model->siesta.grid, grid, 4*sizeof(gint)); /* CURRENT - test for a scaling factor (ie au output instead of angstrom) */ matmat(model->latmat, icell); if (icell[0] > 0.2) { #if DEBUG_SIESTA_DENSITY printf("density grid in au, scale = %f\n", icell[0]); #endif model->siesta.eden_scale = icell[0]; } /* CURRENT - read function data */ switch (type) { case 0: model->siesta.eden = g_malloc(n*grid[0]*sizeof(gdouble)); ptr = model->siesta.eden; break; case 1: model->siesta.epot = g_malloc(n*grid[0]*sizeof(gdouble)); ptr = model->siesta.epot; break; default: ptr = NULL; } f = g_malloc(grid[0]*sizeof(float)); k=0; for (i=0 ; isiesta.eden+k++) = *(f+j); */ /* DEBUG */ /* if (i == n-1) { for (j=0 ; jsiesta.eden+j)); } */ } g_free(f); /* this shouldn't happen (file moved while reading?) */ if (k != n*grid[0]) { gui_text_show(ERROR, "Fatal error reading binary file.\n"); /* g_free(model->siesta.eden); model->siesta.eden = NULL; */ } } } /***********************************************************************/ /* process lattice constant line to produce a cartesian scaling factor */ /***********************************************************************/ #define DEBUG_SIESTA_SCALE 0 gdouble siesta_scale_parse(const gchar *text) { gint num_tokens; gchar **buff; gdouble scale = 1.0; g_assert(text != NULL); #if DEBUG_SIESTA_SCALE printf("input line: %s", text); #endif /* process scale */ buff = tokenize(text, &num_tokens); if (num_tokens > 1) scale *= str_to_float(*(buff+1)); /* process units */ if (num_tokens > 2) if ((g_ascii_strncasecmp(*(buff+2), "au", 2) == 0) || (g_ascii_strncasecmp(*(buff+2), "Bohr", 4) == 0)) scale *= AU2ANG; g_strfreev(buff); #if DEBUG_SIESTA_SCALE printf("output scale: %f\n", scale); #endif return(scale); } /********************************/ /* process a species block line */ /********************************/ void seista_parse_species(const gchar *line, GSList **list) { gint num_tokens; gchar **buff; struct species_pak *species; g_assert(line != NULL); buff = tokenize(line, &num_tokens); if (num_tokens > 2) { species = g_malloc(sizeof(struct species_pak)); species->number = str_to_float(*(buff+1)); species->label = g_strdup(*(buff+2)); *list = g_slist_append(*list, species); } g_strfreev(buff); } /************************************/ /* process a lattice parameter line */ /************************************/ void seista_parse_cell(const gchar *line, gdouble scale, struct model_pak *model) { gint num_tokens; gchar **buff; g_assert(line != NULL); buff = tokenize(line, &num_tokens); if (num_tokens > 5) { model->pbc[0] = scale * str_to_float(*(buff+0)); model->pbc[1] = scale * str_to_float(*(buff+1)); model->pbc[2] = scale * str_to_float(*(buff+2)); model->pbc[3] = D2R*str_to_float(*(buff+3)); model->pbc[4] = D2R*str_to_float(*(buff+4)); model->pbc[5] = D2R*str_to_float(*(buff+5)); } model->construct_pbc = FALSE; /* hack for determining if it's a (GDIS created) surface */ if (fabs(model->pbc[2]) < FRACTION_TOLERANCE) model->periodic = 2; else model->periodic = 3; g_strfreev(buff); } /*********************************/ /* process a lattice vector line */ /*********************************/ void seista_parse_lattice(const gchar *line, gint n, gdouble scale, struct model_pak *model) { gint num_tokens; gchar **buff; g_assert(line != NULL); buff = tokenize(line, &num_tokens); if (num_tokens > 2) { model->latmat[0+n] = scale * str_to_float(*(buff+0)); model->latmat[3+n] = scale * str_to_float(*(buff+1)); model->latmat[6+n] = scale * str_to_float(*(buff+2)); } model->construct_pbc = TRUE; model->periodic = 3; g_strfreev(buff); } /***********************************/ /* process a coordinate block line */ /***********************************/ void seista_parse_coords(const gchar *line, GSList *species_list, gdouble scale, struct model_pak *model) { gint i, num_tokens; gchar **buff; struct species_pak *species; struct core_pak *core; g_assert(line != NULL); buff = tokenize(line, &num_tokens); if (num_tokens > 3) { /* find corresponding number of this element in the species block */ /* NB: list counts from 0, fdf counts from 1 (hence the minus 1) */ if (species_list) { i = abs((gint) str_to_float(*(buff+3)) - 1); species = g_slist_nth_data(species_list, i); core = new_core(species->label, model); /* translucent ghost atom */ if (species->number < 0) { core->ghost = TRUE; core->colour[3] = 0.5; } } else core = new_core("X", model); model->cores = g_slist_prepend(model->cores, core); core->x[0] = str_to_float(*(buff+0)); core->x[1] = str_to_float(*(buff+1)); core->x[2] = str_to_float(*(buff+2)); VEC3MUL(core->x, scale); } g_strfreev(buff); } /***************************/ /* zmatrix line processing */ /***************************/ void seista_parse_zmatrix(const gchar *text, GSList *species_list, gdouble scale, struct model_pak *model) { gint i, num_tokens; gchar **buff; struct core_pak *core; struct species_pak *species_data; static gint state=0; /* keyword processing within zmatrix block */ if (g_ascii_strncasecmp("molecule", text, 8) == 0) { state = 1; if (strstr(text, "frac")) zmat_fractional_set(model->zmatrix); else zmat_cartesian_set(model->zmatrix); return; } if (g_ascii_strncasecmp("variable", text, 8) == 0) { state = 2; return; } if (g_ascii_strncasecmp("constant", text, 8) == 0) { state = 3; return; } if (g_ascii_strncasecmp("fractional", text, 10) == 0) { state = 4; model->fractional = TRUE; return; } if (g_ascii_strncasecmp("cartesian", text, 9) == 0) { state = 4; model->fractional = FALSE; return; } /* data processing within zmatrix block */ switch (state) { case 1: /* process zmatrix coords */ zmat_core_add(text, model->zmatrix); break; case 2: /* process zmatrix variables */ zmat_var_add(text, model->zmatrix); break; case 3: /* process zmatrix constants */ zmat_const_add(text, model->zmatrix); break; case 4: /* fractional input coords */ buff = tokenize(text, &num_tokens); if (num_tokens > 3) { /* find corresponding number of this element in the species block */ /* NB: list counts from 0, fdf counts from 1 (hence the minus 1) */ if (species_list) { i = abs((gint) str_to_float(*(buff)) - 1); species_data = g_slist_nth_data(species_list, i); core = new_core(species_data->label, model); /* translucent ghost atom */ if (species_data->number < 0) { core->ghost = TRUE; core->colour[3] = 0.5; } } else core = new_core("X", model); model->cores = g_slist_prepend(model->cores, core); core->x[0] = scale * str_to_float(*(buff+1)); core->x[1] = scale * str_to_float(*(buff+2)); core->x[2] = scale * str_to_float(*(buff+3)); } g_strfreev(buff); break; } } /******************************/ /* main fdf data block reader */ /******************************/ #define DEBUG_READ_FDF_BLOCK 0 gint read_fdf_block(FILE *fp, struct model_pak *model, gboolean read_inital_cooords) { gint count, num_tokens, block_type; gchar **buff, *line; gdouble scale_lattice = 1.0, scale_coords = 1.0; GSList *species_list=NULL; g_assert(fp != NULL); g_assert(model != NULL); for (;;) { line = file_read_line(fp); /* terminate if NULL returned */ if (!line) break; buff = tokenize(line, &num_tokens); /* get next line if blank */ if (!buff) continue; /* NEW - block/enblock processing */ if (g_ascii_strncasecmp("\%block", line, 6) == 0) { block_type = SIESTA_DEFAULT; if (num_tokens > 1) { if (g_ascii_strncasecmp("ChemicalSpecies", *(buff+1), 15) == 0) block_type = SIESTA_SPECIES; if (g_ascii_strncasecmp("LatticeParameters", *(buff+1), 17) == 0) block_type = SIESTA_LATTICE_PARAMETERS; if (g_ascii_strncasecmp("LatticeVectors", *(buff+1), 14) == 0) block_type = SIESTA_LATTICE_VECTORS; if (g_ascii_strncasecmp("AtomicCoordinates", *(buff+1), 17) == 0) block_type = SIESTA_COORDS; if (g_ascii_strncasecmp("zmatrix", *(buff+1), 7) == 0) { block_type = SIESTA_ZMATRIX; } } #if DEBUG_READ_FDF_BLOCK printf("processing block [type %d]\n", block_type); #endif count = 0; for (;;) { line = file_read_line(fp); if (!line) goto siesta_done_file; if (g_ascii_strncasecmp("\%endblock", line, 9) == 0) { #if DEBUG_READ_FDF_BLOCK printf("end of block.\n"); #endif goto siesta_done_line; } switch (block_type) { case SIESTA_SPECIES: seista_parse_species(line, &species_list); break; case SIESTA_LATTICE_PARAMETERS: seista_parse_cell(line, scale_lattice, model); break; case SIESTA_LATTICE_VECTORS: seista_parse_lattice(line, count, scale_lattice, model); break; case SIESTA_COORDS: seista_parse_coords(line, species_list, scale_coords, model); break; case SIESTA_ZMATRIX: seista_parse_zmatrix(line, species_list, scale_coords, model); break; default: error_table_entry("Unrecognized block encountered.\n"); break; } count++; } /* done block processing - get next line */ continue; } /* cartesian/fractional */ if (read_inital_cooords) { if (g_ascii_strncasecmp("AtomicCoordinatesFormat", line, 23) == 0) { if (g_strrstr(line, "ractional")) model->fractional = TRUE; else { model->fractional = FALSE; /* NEW */ if (g_strrstr(line, "Bohr")) scale_coords = AU2ANG; } } } /* scaling */ if (g_ascii_strncasecmp("LatticeConstant", *buff, 15) == 0) scale_lattice = siesta_scale_parse(line); if (g_ascii_strncasecmp("NumberOfAtoms", *buff, 13) == 0) { if (num_tokens > 1) model->siesta.num_atoms = (gint) str_to_float(*(buff+1)); } if (g_ascii_strncasecmp("SystemLabel", *buff, 11) == 0) { if (num_tokens > 1) { g_free(model->basename); model->basename = g_strdup(*(buff+1)); } } //Terry Siesta options - arrrgh! if (g_ascii_strncasecmp("kgrid_cutoff", *buff, 12) == 0) { if (num_tokens >= 3) //kgrid_cutoff real length_type #othercrap { model->siesta.kgrid_cutoff = str_to_float(*(buff+1)); if (set_length_units(*(buff+2), &model->siesta.kgrid_cutoff_units) !=0) error_table_entry("Units missing kgrid_cutoff.\n"); } else error_table_entry("kgrid_cutoff line corrupt.\n"); } if (g_ascii_strncasecmp("WriteCoorInitial", *buff, 16) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_coor_init = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_coor_init = TRUE; else model->siesta.file_output_write_coor_init = FALSE; } else { error_table_entry("WriteCoorInitial line corrupt\n"); } } if (g_ascii_strncasecmp("WriteCoorStep", *buff, 13) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_coor_step = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_coor_step = TRUE; else model->siesta.file_output_write_coor_step = FALSE; } else error_table_entry("WriteCoorStep line corrupt\n"); } if (g_ascii_strncasecmp("WriteForces", *buff, 11) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_forces = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_forces = TRUE; else model->siesta.file_output_write_forces = FALSE; } else error_table_entry("WriteForcesStep line corrupt\n"); } if (g_ascii_strncasecmp("WriteKpoints", *buff, 12) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_kpoints = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_kpoints = TRUE; else model->siesta.file_output_write_kpoints = FALSE; } else error_table_entry("WriteKpoints line corrupt\n"); } if (g_ascii_strncasecmp("WriteEigenvalues", *buff, 16) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_eigenvalues = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_eigenvalues = TRUE; else model->siesta.file_output_write_eigenvalues = FALSE; } else error_table_entry("WriteEigenvalues line corrupt\n"); } if (g_ascii_strncasecmp("WriteKbands", *buff, 11) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_kbands = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_kbands = TRUE; else model->siesta.file_output_write_kbands = FALSE; } else error_table_entry("WriteKbands line corrupt\n"); } if (g_ascii_strncasecmp("WriteBands", *buff, 10) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_bands = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_bands = TRUE; else model->siesta.file_output_write_bands = FALSE; } else error_table_entry("WriteBands line corrupt\n"); } if (g_ascii_strncasecmp("WriteWaveFunctions", *buff, 18) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_wavefunctions = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_wavefunctions = TRUE; else model->siesta.file_output_write_wavefunctions = FALSE; } else error_table_entry("WriteWaveFunctions line corrupt\n"); } if (g_ascii_strncasecmp("WriteDM", *buff, 7) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_dm = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_dm = TRUE; else model->siesta.file_output_write_dm = FALSE; } else error_table_entry("WriteDM line corrupt\n"); } if (g_ascii_strncasecmp("WriteCoorXmol", *buff, 13) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_coor_xmol = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_coor_xmol = TRUE; else model->siesta.file_output_write_coor_xmol = FALSE; } else error_table_entry("WriteCoorXmol line corrupt\n"); } if (g_ascii_strncasecmp("WriteCoorCerius", *buff, 15) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_coor_cerius = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_coor_cerius = TRUE; else model->siesta.file_output_write_coor_cerius = FALSE; } else error_table_entry("WriteCoorCerius line corrupt\n"); } if (g_ascii_strncasecmp("WriteMDXmol", *buff, 11) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_md_xmol = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_md_xmol = TRUE; else model->siesta.file_output_write_md_xmol = FALSE; } else error_table_entry("WriteMDXmol line corrupt\n"); } if (g_ascii_strncasecmp("WriteMDhistory", *buff, 14) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.file_output_write_md_history = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.file_output_write_md_history = TRUE; else model->siesta.file_output_write_md_history = FALSE; } else error_table_entry("WriteMDhistory line corrupt\n"); } if (g_ascii_strncasecmp("LongOutput", *buff, 10) == 0) { if (num_tokens >= 2) { if ((g_ascii_strncasecmp("true", *buff, 4) == 0) || (g_ascii_strncasecmp(".true", *buff, 5) == 0) ) { model->siesta.long_output = TRUE; model->siesta.file_output_write_coor_step = TRUE; model->siesta.file_output_write_forces = TRUE; model->siesta.file_output_write_kpoints = TRUE; model->siesta.file_output_write_eigenvalues = TRUE; model->siesta.file_output_write_bands = TRUE; model->siesta.file_output_write_kbands = TRUE; model->siesta.file_output_write_wavefunctions = TRUE; model->siesta.file_output_write_mullikenpop = 1.0; } else { model->siesta.long_output = FALSE; model->siesta.file_output_write_coor_step = FALSE; model->siesta.file_output_write_forces = FALSE; model->siesta.file_output_write_kpoints = FALSE; model->siesta.file_output_write_eigenvalues = FALSE; model->siesta.file_output_write_bands = FALSE; model->siesta.file_output_write_kbands = FALSE; model->siesta.file_output_write_wavefunctions = FALSE; model->siesta.file_output_write_mullikenpop = 0.0; } } else { model->siesta.long_output = FALSE; model->siesta.file_output_write_coor_step = FALSE; model->siesta.file_output_write_forces = FALSE; model->siesta.file_output_write_kpoints = FALSE; model->siesta.file_output_write_eigenvalues = FALSE; model->siesta.file_output_write_bands = FALSE; model->siesta.file_output_write_kbands = FALSE; model->siesta.file_output_write_wavefunctions = FALSE; model->siesta.file_output_write_mullikenpop = 0.0; error_table_entry("LongOutput line corrupt\n"); } } if (g_ascii_strncasecmp("PAO.SplitNorm", *buff, 13) == 0) { if (num_tokens >= 2) model->siesta.split_zeta_norm = str_to_float(*(buff+1)); else error_table_entry("PAO.SplitNorm line corrupt : using defaults.\n"); } if (g_ascii_strncasecmp("PAO.BasisSize", *buff, 13) == 0) { if (num_tokens == 2) { if (g_ascii_strncasecmp("SZ", *(buff+1), 2) == 0) model->siesta.basis_set = SZ_ZETA; if (g_ascii_strncasecmp("SZP", *(buff+1), 3) == 0) model->siesta.basis_set = SZP_ZETA; if (g_ascii_strncasecmp("DZ", *(buff+1), 2) == 0) model->siesta.basis_set = DZ_ZETA; if (g_ascii_strncasecmp("DZP", *(buff+1), 3) == 0) model->siesta.basis_set = DZP_ZETA; } else { model->siesta.basis_set = DZP_ZETA; error_table_entry("PAO.BasisSize line corrupt : using defaults\n"); } } if (g_ascii_strncasecmp("PAO.EnergyShift", *buff, 15) == 0) { if (num_tokens >= 3) { model->siesta.energy_shift = str_to_float(*(buff+1)); if (set_energy_units(*(buff+2), &model->siesta.energy_shift_units) != 0) error_table_entry("PAO.EnergyShift line corrupt : using defaults\n"); } else error_table_entry("PAO.EnergyShift line corrupt : using defaults\n"); } if (g_ascii_strncasecmp("SolutionMethod", *buff, 14) == 0) { if (g_ascii_strncasecmp("OrderN", *(buff+1), 6) == 0) model->siesta.solution_method_type = ORDERN_SOLMETH; else if (g_ascii_strncasecmp("diagon", *(buff+1), 6) == 0) model->siesta.solution_method_type = DIAGON_SOLMETH; else error_table_entry("SolutionMethod line corrupt\n"); } if (g_ascii_strncasecmp("XC.functional", *buff, 13) == 0) { if (g_ascii_strncasecmp("LDA", *(buff+1), 3) == 0) model->siesta.xc_functional_type = LDA_XCFUNC; else if (g_ascii_strncasecmp("LSD", *(buff+1), 3) == 0) model->siesta.xc_functional_type = LSD_XCFUNC; else if (g_ascii_strncasecmp("GGA", *(buff+1), 3) == 0) model->siesta.xc_functional_type = GGA_XCFUNC; else error_table_entry("XC.functional line corrupt\n"); } if (g_ascii_strncasecmp("XC.authors", *buff, 10) == 0) { if (g_ascii_strncasecmp("PZ", *(buff+1), 2) == 0) model->siesta.xc_author_type = PZ_XCAUTH; else if (g_ascii_strncasecmp("CA", *(buff+1), 2) == 0) model->siesta.xc_author_type = CA_XCAUTH; else if (g_ascii_strncasecmp("PW92", *(buff+1), 4) == 0) model->siesta.xc_author_type = PW_XCAUTH; else if (g_ascii_strncasecmp("PBE", *(buff+1), 3) == 0) model->siesta.xc_author_type = PBE_XCAUTH; else if (g_ascii_strncasecmp("RPBE", *(buff+1), 4) == 0) model->siesta.xc_author_type = RPBE_XCAUTH; else if (g_ascii_strncasecmp("LYP", *(buff+1), 4) == 0) model->siesta.xc_author_type = LYP_XCAUTH; else error_table_entry("XC.authors line corrupt\n"); } //mesh_cutoff -> MeshCutoff real siesta_unit_type if (g_ascii_strncasecmp("MeshCutoff", *buff, 10) == 0) { if (num_tokens >= 3) { model->siesta.mesh_cutoff = str_to_float(*(buff+1)); if (set_energy_units(*(buff+2), &model->siesta.mesh_cutoff_units) != 0) error_table_entry("MeshCutoff UNITS corrupt : using default\n"); } else error_table_entry("MeshCutoff line corrupt : using defaults\n"); } if (g_ascii_strncasecmp("Harris_functional", *buff, 17) == 0) { if (num_tokens >=2) //Harris_functional boolean #crap { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.harris_functional = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.harris_functional = TRUE; else model->siesta.harris_functional = FALSE; } else error_table_entry("Harris_functional line corrupt\n"); } if (g_ascii_strncasecmp("ElectronicTemperature", *buff, 21) == 0) { if (num_tokens > 2) //temp real units #othercrap { model->siesta.electronic_temperature = str_to_float(*(buff+1)); if (set_temperature_units(*(buff+2), &model->siesta.electronic_temperature_units) != 0) error_table_entry("ElectronicTemperature UNITS corrupt : using default\n"); } else if (num_tokens == 2) { model->siesta.electronic_temperature = str_to_float(*(buff+1)); //default temp units? KELVIN FOR NOW! model->siesta.electronic_temperature_units = K_TEMP; error_table_entry("Units missing ElectronicTemperature : using Kelvin\n"); } else error_table_entry("ElectronicTemperature line corrupt : using defaults\n"); } //spin_polarised -> SpinPolarized true||.true. if (g_ascii_strncasecmp("SpinPolarized", *buff, 13) == 0) { if (num_tokens >= 2) //SpinPolarized boolean #othercrap { if (g_ascii_strncasecmp("true", *(buff+1), 4) == 0) model->siesta.spin_polarised = TRUE; else if (g_ascii_strncasecmp(".true", *(buff+1), 5) == 0) model->siesta.spin_polarised = TRUE; else model->siesta.spin_polarised = FALSE; } else error_table_entry("SpinPolarized line corrupt : using defaults\n"); } //SIESTA - SCF OPTIONS // -----> right no_of_cycles -> MaxSCFIterations if (g_ascii_strncasecmp("MaxSCFIterations", *buff, 16) == 0) { if (num_tokens >= 2) model->siesta.no_of_cycles = str_to_float(*(buff+1)); else error_table_entry("MaxSCFIterations line corrupt : using defaults\n"); } if (g_ascii_strncasecmp("DM.NumberPulay", *buff, 14) == 0) { if (num_tokens >= 2) model->siesta.no_of_pulay_matrices = str_to_float(*(buff+1)); else error_table_entry("DM.NumberPulay line corrupt : using defaults\n"); } //mixing_weight -> DM.MixingWeight float #comments if (g_ascii_strncasecmp("DM.MixingWeight", *buff, 15) == 0) { if (num_tokens >= 2) model->siesta.mixing_weight = str_to_float(*(buff+1)); else error_table_entry("DM.MixingWeight line corrupt : using defaults\n"); } //run_type - enum { SINGLE_POINT, OPTIMISATION, MOLECULAR_DYNAMICS, PHONON_CALCULATION }; if (g_ascii_strncasecmp("MD.TypeOfRun", *buff, 12) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("CG", *(buff+1), 2) == 0) { //soo is it a single point or an optimisation? //have to check for number of steps if (model->siesta.number_of_steps == 0) model->siesta.run_type = SINGLE_POINT; else model->siesta.run_type = OPTIMISATION; } else if (g_ascii_strncasecmp("Verlet", *(buff+1), 6) == 0) { model->siesta.run_type = MOLECULAR_DYNAMICS; model->siesta.md_type_of_run = VERLET_MDRUN; } else if (g_ascii_strncasecmp("Nose", *(buff+1), 4) == 0) { model->siesta.run_type = MOLECULAR_DYNAMICS; model->siesta.md_type_of_run = NOSE_MDRUN; } else if (g_ascii_strncasecmp("ParrinelloRahman", *(buff+1), 16) == 0) { model->siesta.run_type = MOLECULAR_DYNAMICS; model->siesta.md_type_of_run = PARRINELLOPAHMAN_MDRUN; } else if (g_ascii_strncasecmp("NoseParrinelloRahman", *(buff+1), 20) == 0) { model->siesta.run_type = MOLECULAR_DYNAMICS; model->siesta.md_type_of_run = NOSEPARRINELLOPAHMAN_MDRUN; } else if (g_ascii_strncasecmp("Anneal", *(buff+1), 6) == 0) { model->siesta.run_type = MOLECULAR_DYNAMICS; model->siesta.md_type_of_run = ANNEAL_MDRUN; } else if (g_ascii_strncasecmp("FC", *(buff+1), 2) == 0) { model->siesta.run_type = PHONON_CALCULATION; model->siesta.md_type_of_run = FC_MDRUN; } else if (g_ascii_strncasecmp("Phonon", *(buff+1), 6) == 0) { model->siesta.run_type = PHONON_CALCULATION; model->siesta.md_type_of_run = PHONON_MDRUN; } } else error_table_entry("MD.TypeOfRun line corrupt"); } //model->siesta.md_variable_cell if (g_ascii_strncasecmp("MD.VariableCell", *buff, 15) == 0) { if (num_tokens >=2) { if (set_true_false(*(buff+1), &model->siesta.md_variable_cell) != 0) error_table_entry("Failed to set MD.VariableCell\n"); } else error_table_entry("MD.VariableCell line corrupt : using defaults\n"); } //model->siesta.md_num_cg_steps if (g_ascii_strncasecmp("MD.NumCGsteps", *buff, 13) == 0) { if (num_tokens >= 2) model->siesta.md_num_cg_steps = str_to_float(*(buff+1)); else error_table_entry("MD.NumCGsteps line corrupt\n"); } //model->siesta.md_max_cg_displacement if (g_ascii_strncasecmp("MD.MaxCGDispl", *buff, 13) == 0) { if (num_tokens >= 3) { model->siesta.md_max_cg_displacement = str_to_float(*(buff+1)); if (set_length_units(*(buff+2), &model->siesta.md_max_cg_displacement_units) != 0) error_table_entry("MD.MaxCGDispl units unable to be set\n"); } else error_table_entry("MD.MaxCGDispl line corrupt\n"); } //model->siesta.md_precondition_variable_cell if (g_ascii_strncasecmp("MD.PreconditionVariableCell", *buff, 27) == 0) { if (num_tokens >= 3) { model->siesta.md_precondition_variable_cell = str_to_float(*(buff+1)); if (set_length_units(*(buff+2), &model->siesta.md_precondition_variable_cell_units) != 0) error_table_entry("MD.PreconditionVariableCell units unable to be set\n"); } else error_table_entry("MD.PreconditionVariableCell line corrupt\n"); } //model->siesta.md_max_force_tol if (g_ascii_strncasecmp("MD.MaxForceTol", *buff, 14) == 0) { if (num_tokens >= 3) { model->siesta.md_max_force_tol = str_to_float(*(buff+1)); if (set_force_units(*(buff+2), &model->siesta.md_max_force_tol_units) != 0) error_table_entry("MD.MaxForceTol units unable to be set\n"); } else error_table_entry("MD.MaxForceTol line corrupt\n"); } //model->siesta.md_max_stress_tol if (g_ascii_strncasecmp("MD.MaxStressTol", *buff, 14) == 0) { if (num_tokens >= 3) { model->siesta.md_max_stress_tol = str_to_float(*(buff+1)); if (set_pressure_units(*(buff+2), &model->siesta.md_max_stress_tol_units) != 0) error_table_entry("MD.MaxStressTol units unable to be set\n"); } else error_table_entry("MD.MaxStressTol line corrupt\n"); } //model->siesta.md_inital_time_step if (g_ascii_strncasecmp("MD.InitialTimeStep", *buff, 18) == 0) { if (num_tokens >= 2) model->siesta.md_inital_time_step = str_to_float(*(buff+1)); else error_table_entry("MD.InitialTimeStep line corrupt\n"); } //model->siesta.md_final_time_step if (g_ascii_strncasecmp("MD.FinalTimeStep", *buff, 16) == 0) { if (num_tokens >= 2) model->siesta.md_final_time_step = str_to_float(*(buff+1)); else error_table_entry("MD.FinalTimeStep line corrupt\n"); } //model->siesta.timestep if (g_ascii_strncasecmp("MD.LengthTimeStep", *buff, 17) == 0) { if (num_tokens >= 3) { model->siesta.timestep = str_to_float(*(buff+1)); if (set_time_units(*(buff+2), &model->siesta.timestep_units) != 0) error_table_entry("*ERROR* - MD.LengthTimeStep units corrupt\n"); } else if (num_tokens == 2) { error_table_entry("MD.LengthTimeStep units MISSING\n"); model->siesta.timestep = str_to_float(*(buff+1)); set_time_units("", &model->siesta.timestep_units); } } //model->siesta.md_inital_temperature if (g_ascii_strncasecmp("MD.InitialTemperature", *buff, 21) == 0) { if (num_tokens >= 3) { model->siesta.md_inital_temperature = str_to_float(*(buff+1)); if (set_temperature_units(*(buff+2), &model->siesta.md_inital_temperature_units) != 0) error_table_entry("MD.InitialTemperature units corrupt\n"); } else error_table_entry("MD.InitialTemperature line corrupt\n"); } //model->siesta.md_target_temperature if (g_ascii_strncasecmp("MD.TargetTemperature", *buff, 20) == 0) { if (num_tokens >= 3) { model->siesta.md_target_temperature = str_to_float(*(buff+1)); if (set_temperature_units(*(buff+2), &model->siesta.md_target_temperature_units) != 0) error_table_entry("MD.TargetTemperature units corrupt\n"); } else error_table_entry("MD.TargetTemperature line corrupt\n"); } //model->siesta.md_quench if (g_ascii_strncasecmp("MD.Quench", *buff, 9) == 0) { if (num_tokens >= 2) { if (set_true_false(*(buff+1), &model->siesta.md_quench) != 0) error_table_entry("MD.Quench boolean corrupt\n"); } else error_table_entry("MD.Quench line corrupt\n"); } if (g_ascii_strncasecmp("MD.TargetPressure", *buff, 20) == 0) { if (num_tokens >= 3) { model->siesta.md_target_pressure = str_to_float(*(buff+1)); if (set_pressure_units(*(buff+2), &model->siesta.md_target_pressure_units) != 0) error_table_entry("MD.TargetPressure units corrupt\n"); } else error_table_entry("MD.TargetPressure line corrupt\n"); } #if OLD_STYLE if (g_ascii_strncasecmp("\%block", *buff, 6) == 0 && num_tokens > 1) { if (g_ascii_strncasecmp("MD.TargetStress", *(buff+1), 15) == 0) { for (;;) { g_strfreev(buff); buff = get_tokenized_line(fp, &num_tokens); if (!buff) break; if (g_ascii_strncasecmp("\%endblock", *buff, 9) == 0) break; if (num_tokens == 6) { model->siesta.md_target_stress_xx = str_to_float(*(buff)); model->siesta.md_target_stress_yy = str_to_float(*(buff+1)); model->siesta.md_target_stress_zz = str_to_float(*(buff+2)); model->siesta.md_target_stress_xy = str_to_float(*(buff+3)); model->siesta.md_target_stress_xz = str_to_float(*(buff+4)); model->siesta.md_target_stress_yz = str_to_float(*(buff+5)); } else { //invalid stress block text = error_table_entry("*ERROR* - MD.TargetStress block corrupt\n"); gui_text_show(ERROR, text); g_free(text); } } } } #endif if (g_ascii_strncasecmp("MD.NoseMass", *buff, 11) == 0) { if (num_tokens >= 3) { model->siesta.md_nose_mass = str_to_float(*(buff+1)); if (set_mominert_units(*(buff+2), &model->siesta.md_nose_mass_units) != 0) error_table_entry("*ERROR* - MD.NoseMass units corrupt\n"); } else error_table_entry("*ERROR* - MD.NoseMass line corrupt\n"); } if (g_ascii_strncasecmp("MD.ParrinelloRahmanMass", *buff, 23) == 0) { if (num_tokens >= 3) { model->siesta.md_parrinello_rahman_mass = str_to_float(*(buff+1)); if (set_mominert_units(*(buff+2), &model->siesta.md_parrinello_rahman_mass_units) != 0) error_table_entry("MD.ParrinelloRahmanMass units corrupt\n"); } else error_table_entry("MD.ParrinelloRahmanMass line corrupt\n"); } if (g_ascii_strncasecmp("MD.AnnealOption", *buff, 23) == 0) { if (num_tokens >= 2) { if (g_ascii_strncasecmp("Temperature", *(buff+1), 11) == 0) model->siesta.md_anneal_option = TEMPERATURE_ANNEAL; else if (g_ascii_strncasecmp("Pressure", *(buff+1), 8) == 0) model->siesta.md_anneal_option = PRESSURE_ANNEAL; else if (g_ascii_strncasecmp("TemperatureandPressure", *(buff+1), 22) == 0) model->siesta.md_anneal_option = TEMP_AND_PRESS_ANNEAL; else error_table_entry("MD.AnnealOption string corrupt\n"); } else error_table_entry("MD.AnnealOption line corrupt\n"); } if (g_ascii_strncasecmp("MD.TauRelax", *buff, 11) == 0) { if (num_tokens >= 3) { model->siesta.md_tau_relax = str_to_float(*(buff+1)); if (set_time_units(*(buff+2), &model->siesta.md_tau_relax_units) != 0) error_table_entry("MD.TauRelax units corrupt\n"); } else error_table_entry("MD.TauRelax line corrupt\n"); } if (g_ascii_strncasecmp("MD.BulkModulus", *buff, 14) == 0) { if (num_tokens >= 3) { model->siesta.md_bulk_modulus = str_to_float(*(buff+1)); if (set_pressure_units(*(buff+2), &model->siesta.md_bulk_modulus_units) != 0) error_table_entry("MD.BulkModulus units corrupt\n"); } else error_table_entry("MD.BulkModulus line corrupt\n"); } if (g_ascii_strncasecmp("MD.FCDispl", *buff, 10) == 0) { if (num_tokens >= 3) { model->siesta.md_fc_displ = str_to_float(*(buff+1)); if (set_length_units(*(buff+2), &model->siesta.md_fc_displ_units) != 0) error_table_entry("MD.FCDispl units corrupt\n"); } else error_table_entry("MD.FCDispl line corrupt\n"); } if (g_ascii_strncasecmp("MD.FCfirst", *buff, 10) == 0) { if (num_tokens >= 2) model->siesta.md_fc_first = str_to_float(*(buff+1)); else error_table_entry("MD.FCfirst line corrupt\n"); } if (g_ascii_strncasecmp("MD.FClast", *buff, 9) == 0) { if (num_tokens >= 2) model->siesta.md_fc_last = str_to_float(*(buff+1)); else error_table_entry("MD.FClast line corrupt\n"); } /* zmatrix */ if (g_ascii_strncasecmp("zm.unitslength", *buff, 14) == 0) { if (num_tokens) { if (g_ascii_strncasecmp("ang", *(buff+1), 3) == 0) zmat_distance_units_set(model->zmatrix, ANGSTROM); else zmat_distance_units_set(model->zmatrix, BOHR); } } if (g_ascii_strncasecmp("zm.unitsangle", *buff, 13) == 0) { if (num_tokens) { if (g_ascii_strncasecmp("deg", *(buff+1), 3) == 0) zmat_angle_units_set(model->zmatrix, DEGREES); else zmat_angle_units_set(model->zmatrix, RADIANS); } } #if OLD_STYLE if (read_inital_cooords) { if (g_ascii_strncasecmp("AtomicCoordinates", *(buff+1), 17) == 0) { } } #endif // exit condition - "************************** End of input data file" if (g_ascii_strncasecmp("**************************", *buff, 26) == 0) { if (g_ascii_strncasecmp("End", *(buff+1), 3) == 0) { //leave routine -> break the for loop g_strfreev(buff); free_slist(species_list); return 1; } } /* loop cleanup */ siesta_done_line:; g_strfreev(buff); g_free(line); } siesta_done_file:; /* NEW - process zmatrix cores */ zmat_type(model->zmatrix, species_list); zmat_process(model->zmatrix, model); /* free the species list */ free_slist(species_list); return 0; } /**************************/ /* fdf input file reading */ /**************************/ #define DEBUG_READ_FDF 0 gint read_fdf(gchar *filename, struct model_pak *model) { gint i; gchar *text; FILE *fp; /* checks */ g_return_val_if_fail(model != NULL, 1); g_return_val_if_fail(filename != NULL, 2); fp = fopen(filename, "rt"); if (!fp) return(3); error_table_clear(); /* terry's mod */ read_fdf_block(fp, model, TRUE); /* check cores */ model->cores = g_slist_reverse(model->cores); i = g_slist_length(model->cores); if (model->siesta.num_atoms != i) { text = g_strdup_printf("Inconsistancy reading %s: expected %d cores, but found %d.\n", filename, model->siesta.num_atoms, i); gui_text_show(ERROR, text); g_free(text); } /* model setup */ strcpy(model->filename, filename); model_prep(model); /* NB: this relies on having a valid latmat */ seista_density_files(filename, model); error_table_print_all(); return(0); } /*******************************************/ /* read single SIESTA output configuration */ /*******************************************/ #define DEBUG_READ_SOUT 0 gint read_sout_block(FILE *fp, struct model_pak *model) { gint i, num_tokens; gchar **temp, **buff, line[LINELEN]; GString *title, *energy_string, *grad_string, *pressure_string; GSList *clist; struct core_pak *core; clist = model->cores; /* get 1st line of coords */ if (fgetline(fp, line)) return(1); buff = tokenize(line, &num_tokens); while (num_tokens > 4) { if (clist) { core = 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] = str_to_float(*(buff+0)); core->x[1] = str_to_float(*(buff+1)); core->x[2] = str_to_float(*(buff+2)); /* get next line */ g_strfreev(buff); if (fgetline(fp, line)) return(2); buff = tokenize(line, &num_tokens); } g_strfreev(buff); /* attempt to get the lattice matrix */ /* if not found, use the initial values */ /* removed the following line as if initial structure used LatticeVectors, this stopped the pbcs being constructed */ /* model->construct_pbc = FALSE;*/ while (!fgetline(fp, line)) { /* if next frame encoutered - assume there is no outcell data & exit */ if (g_ascii_strncasecmp(line, "siesta: Begin CG move", 28) == 0) break; /* new version of siesta doesn't have the siesta: */ if (g_ascii_strncasecmp(line, " Begin CG move", 41) == 0) break; /* Read in MD steps from siesta: */ if (g_ascii_strncasecmp(line, " Begin MD step", 41) == 0) break; /* acquire cell lengths */ if (g_ascii_strncasecmp(line, "outcell: Cell vector modules", 28) == 0) { model->construct_pbc = FALSE; temp = g_strsplit(line, ":", 3); buff = tokenize(*(temp+2), &num_tokens); g_strfreev(temp); if (num_tokens > 2) { model->pbc[0] = str_to_float(*(buff+0)); model->pbc[1] = str_to_float(*(buff+1)); model->pbc[2] = str_to_float(*(buff+2)); } else printf("Unexpected data format reading cell lengths.\n"); g_strfreev(buff); } /* acquire cell angles */ if (g_ascii_strncasecmp(line, "outcell: Cell angles", 20) == 0) { temp = g_strsplit(line, ":", 3); buff = tokenize(*(temp+2), &num_tokens); g_strfreev(temp); if (num_tokens > 2) { /* get the angles (in radians) */ model->pbc[3] = D2R*str_to_float(*(buff+0)); model->pbc[4] = D2R*str_to_float(*(buff+1)); model->pbc[5] = D2R*str_to_float(*(buff+2)); } else printf("Unexpected data format reading cell angles.\n"); g_strfreev(buff); /* no more data - break out */ /* break; */ } if (g_ascii_strncasecmp(line, "siesta: E_KS(eV)", 16) == 0) { buff = tokenize(line, &num_tokens); model->siesta.energy = str_to_float(*(buff+3)); model->siesta.have_energy = TRUE; g_strfreev(buff); } if (g_strrstr(line, "constrained") != NULL) { buff = tokenize(line, &num_tokens); model->siesta.max_grad = str_to_float(*(buff+1)); model->siesta.have_max_grad = TRUE; g_strfreev(buff); } if (g_ascii_strncasecmp(line, "siesta: Pressure (static)", 25) == 0) { for (i=0; i<4; i++) { if (fgetline(fp, line)) return(2); } buff = tokenize(line, &num_tokens); pressure_string = g_string_new(""); g_string_append_printf(pressure_string, "%.4f %s", str_to_float(*(buff+1)), *(buff+3)); property_add_ranked(5, "Final Pressure", pressure_string->str, model); g_string_free(pressure_string, TRUE); } } g_free(model->title); /* TODO read in units from sot file */ title = g_string_new(""); if (model->siesta.have_energy) { energy_string = g_string_new(""); g_string_append_printf(energy_string, "%.5f eV", model->siesta.energy); property_add_ranked(3, "Energy", energy_string->str, model); g_string_free(energy_string, TRUE); g_string_append_printf(title, "E"); g_string_append_printf(title, " = %.4f eV, ", model->siesta.energy); } if (model->siesta.have_max_grad) { grad_string = g_string_new(""); g_string_append_printf(grad_string, "%.4f eV/A", model->siesta.max_grad); property_add_ranked(4, "Maximum Gradient", grad_string->str, model); g_string_free(grad_string, TRUE); g_string_append_printf(title, "max grad = %.5f", model->siesta.max_grad); } model->title = g_strdup(title->str); g_string_free(title, TRUE); return(0); } /*******************************/ /* SIESTA output frame reading */ /*******************************/ gint read_sout_frame(FILE *fp, struct model_pak *model) { return(read_sout_block(fp, model)); } /********************************/ /* Read in a SIESTA output file */ /********************************/ /* could almost use the read_fdf() code, as SIESTA spits out a copy */ /* however, we do need to know the number of frames for animation... */ gint read_sout(gchar *filename, struct model_pak *model) { gint frame=0; gchar **buff, line[LINELEN], *text; FILE *fp; gint num_tokens, i; gdouble scale = 1.0; fp = fopen(filename, "rt"); if (!fp) return(1); fgetline(fp, line); while (g_ascii_strncasecmp("************************** Dump of input data file", line, 50) != 0) fgetline(fp, line); //should be at the input file params -- dont want inital co-ords. //read_fdf_block(fp,model, TRUE); //should magicly return at the ====== while (!fgetline(fp, line)) { //---------------START TERRY COMMENT DUE TO read_fdf_pak_filler()---------- /* lattice constant scaling */ if (g_ascii_strncasecmp("LatticeConstant", line, 15) == 0) { scale = siesta_scale_parse(line); continue; } /* default cell dimensions */ if (g_ascii_strncasecmp("\%block LatticeParameters", line, 24) == 0) { if (fgetline(fp, line)) return(2); buff = tokenize(line, &num_tokens); if (num_tokens > 5) { model->pbc[0] = scale*str_to_float(*(buff+0)); model->pbc[1] = scale*str_to_float(*(buff+1)); model->pbc[2] = scale*str_to_float(*(buff+2)); model->pbc[3] = D2R*str_to_float(*(buff+3)); model->pbc[4] = D2R*str_to_float(*(buff+4)); model->pbc[5] = D2R*str_to_float(*(buff+5)); model->construct_pbc = FALSE; model->periodic = 3; } g_strfreev(buff); } else if (g_ascii_strncasecmp("\%block LatticeVectors", line, 20) == 0) { for (i=0; i<3; i++) { if (fgetline(fp, line)) return(2); buff = tokenize(line, &num_tokens); if (num_tokens == 3) { model->latmat[0+i] = scale*str_to_float(*(buff+0)); model->latmat[3+i] = scale*str_to_float(*(buff+1)); model->latmat[6+i] = scale*str_to_float(*(buff+2)); } g_strfreev(buff); model->construct_pbc = TRUE; model->periodic = 3; } } //---------------END TERRY COMMENT DUE TO read_fdf_pak_filler()---------- /* functional */ if (g_ascii_strncasecmp(line, "xc_check: Exchange-correlation functional:", 42) == 0) { if (fgetline(fp, line)) return(2); if (g_ascii_strncasecmp(line, "xc_check: GGA Perdew, Burke & Ernzerhof 1996", 44) == 0) property_add_ranked(7, "Functional", "PBE", model); } /* k-points */ if (g_ascii_strncasecmp(line, "siesta: k-grid: Number of k-points", 34) == 0) { buff = g_strsplit(line, "=", 2); g_strchug(g_strchomp(*(buff+1))); property_add_ranked(10, "K-points", *(buff+1), model); g_strfreev(buff); } /* Mesh cutoff */ if (g_ascii_strncasecmp(line, "redata: Mesh Cutoff", 19) == 0) { buff = g_strsplit(line, "=", 2); text = format_value_and_units(*(buff+1), 2); property_add_ranked(6, "Mesh Cutoff", text, model); g_free(text); g_strfreev(buff); } /* Energy Shift */ if (g_ascii_strncasecmp(line, "SPLIT: energy shift", 19) == 0) { buff = g_strsplit(line, "=", 2); text = format_value_and_units(*(buff+1), 6); property_add_ranked(6, "Energy Shift", text, model); g_free(text); g_strfreev(buff); } /* coordinates */ if (model->siesta.md_type_of_run != FC_MDRUN) { if (g_ascii_strncasecmp(line, "outcoor: Atomic coordinates", 27) == 0) { if (g_strrstr(line, "Ang") != NULL) model->fractional = FALSE; else if (g_strrstr(line, "fractional") != NULL) model->fractional = TRUE; else if (g_strrstr(line, "Bohr") != NULL) model->coord_units = BOHR; else { gui_text_show(ERROR, "unexpected coordinate type\n"); return(2); } add_frame_offset(fp, model); read_sout_block(fp, model); frame++; } } /* seems to increment the number of frames but doesn't set the animation pointer! */ /* line only occurs at end of a run that runs out of cycles so commenting out! */ /* if (g_ascii_strncasecmp(line, "outcoor: Final", 14) == 0) frame++;*/ } if (model->siesta.md_type_of_run == FC_MDRUN) { //something went wrong error_table_entry("MD.TypeOfRun == FC -> so inital co-ords from input shown\n"); } /* 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); } gint set_energy_units(gchar * buff, enum siesta_energy_measurement_type * siesta_unit_type) { gint allok = 0; if (g_ascii_strncasecmp("J", buff, 1) == 0) *siesta_unit_type = J_ENRG; else if (g_ascii_strncasecmp("erg", buff, 3) == 0) *siesta_unit_type = ERG_ENRG; else if (g_ascii_strncasecmp("eV", buff, 2) == 0) *siesta_unit_type = EV_ENRG; else if (g_ascii_strncasecmp("meV", buff, 3) == 0) *siesta_unit_type = MEV_ENRG; else if (g_ascii_strncasecmp("Ry", buff, 2) == 0) *siesta_unit_type = RY_ENRG; else if (g_ascii_strncasecmp("mRy", buff, 3) == 0) *siesta_unit_type = MRY_ENRG; else if (g_ascii_strncasecmp("Hartree", buff, 7) == 0) *siesta_unit_type = HARTREE_ENRG; else if (g_ascii_strncasecmp("K", buff, 1) == 0) *siesta_unit_type = KO_ENRG; else if (g_ascii_strncasecmp("kcal/mol", buff, 8) == 0) *siesta_unit_type = KCALMOLE_ENRG; else if (g_ascii_strncasecmp("mHartree", buff, 8) == 0) *siesta_unit_type = MHARTREE_ENRG; else if (g_ascii_strncasecmp("kJ/mol", buff, 6) == 0) *siesta_unit_type = KJMOL_ENRG; else if (g_ascii_strncasecmp("Hz", buff, 3) == 0) *siesta_unit_type = HZ_ENRG; else if (g_ascii_strncasecmp("THz", buff, 3) == 0) *siesta_unit_type = THZ_ENRG; else if (g_ascii_strncasecmp("cm-1", buff, 4) == 0) *siesta_unit_type = CM_ENRG; else if (g_ascii_strncasecmp("cm**-1", buff, 6) == 0) *siesta_unit_type = CMCM_ENRG; else { //no match - set default *siesta_unit_type = RY_ENRG; allok = 1; } return allok; } gint set_temperature_units(gchar * buff, enum siesta_temperature_measurement_type * siesta_unit_type) { gint allok = 0; if (g_ascii_strncasecmp("K", buff, 1) == 0) *siesta_unit_type = K_TEMP; else if (g_ascii_strncasecmp("C", buff, 1) == 0) *siesta_unit_type = C_TEMP; else { //no match - set default *siesta_unit_type = K_TEMP; allok = 1; } return allok; } gint set_mass_units(gchar * buff, enum siesta_temperature_measurement_type * siesta_unit_type) { //KG_MASS, G_MASS, AMU_MASS // gint allok = 0; if (g_ascii_strncasecmp("Kg", buff, 2) == 0) *siesta_unit_type = KG_MASS; else if (g_ascii_strncasecmp("g", buff, 1) == 0) *siesta_unit_type = G_MASS; else if (g_ascii_strncasecmp("amu", buff, 3) == 0) *siesta_unit_type = AMU_MASS; else { //no match - set default *siesta_unit_type = AMU_MASS; allok = 1; } return allok; } gint set_length_units(gchar * buff, enum siesta_length_measurement_type * siesta_unit_type) { //M_LEN, CM_LEN, NM_LEN, ANG_LEN, BOHR_LEN gint allok = 0; if (g_ascii_strncasecmp("m", buff, 1) == 0) *siesta_unit_type = M_LEN; else if (g_ascii_strncasecmp("cm", buff, 2) == 0) *siesta_unit_type = CM_LEN; else if (g_ascii_strncasecmp("nm", buff, 2) == 0) *siesta_unit_type = NM_LEN; else if (g_ascii_strncasecmp("Ang", buff, 3) == 0) *siesta_unit_type = ANG_LEN; else if (g_ascii_strncasecmp("Bohr", buff, 4) == 0) *siesta_unit_type = BOHR_LEN; else { //no match - set default *siesta_unit_type = BOHR_LEN; allok = 1; } return allok; } gint set_time_units(gchar * buff, enum siesta_time_measurement_type * siesta_unit_type) { //S_TIME, FS_TIME, PS_TIME, NS_TIME gint allok = 0; if (g_ascii_strncasecmp("S", buff, 1) == 0) *siesta_unit_type = S_TIME; else if (g_ascii_strncasecmp("fs", buff, 2) == 0) *siesta_unit_type = FS_TIME; else if (g_ascii_strncasecmp("ps", buff, 2) == 0) *siesta_unit_type = PS_TIME; else if (g_ascii_strncasecmp("ns", buff, 2) == 0) *siesta_unit_type = NS_TIME; else { //no match - set default *siesta_unit_type = FS_TIME; allok = 1; } return allok; } gint set_pressure_units(gchar * buff, enum siesta_pressure_measurement_type * siesta_unit_type) { // P_PRES, MPA_PRES, GPA_PRES, ATM_PRES, BAR_PRES, MBAR_PRES, RYBOHR3_PRES, EVANG3_PRES gint allok = 0; if (g_ascii_strncasecmp("Pa", buff, 2) == 0) *siesta_unit_type = P_PRES; else if (g_ascii_strncasecmp("MPa", buff, 3) == 0) *siesta_unit_type = MPA_PRES; else if (g_ascii_strncasecmp("GPa", buff, 3) == 0) *siesta_unit_type = GPA_PRES; else if (g_ascii_strncasecmp("atm", buff, 3) == 0) *siesta_unit_type = ATM_PRES; else if (g_ascii_strncasecmp("bar", buff, 3) == 0) *siesta_unit_type = BAR_PRES; else if (g_ascii_strncasecmp("Mbar", buff, 4) == 0) *siesta_unit_type = MBAR_PRES; else if (g_ascii_strncasecmp("Ry/Bohr**3", buff, 10) == 0) *siesta_unit_type = RYBOHR3_PRES; else if (g_ascii_strncasecmp("eV/Ang**3", buff, 9) == 0) *siesta_unit_type = EVANG3_PRES; else { //no match - set default *siesta_unit_type = P_PRES; allok = 1; } return allok; } gint set_force_units(gchar * buff, enum siesta_force_measurement_type * siesta_unit_type) { //N_FORCE, EVANG_FORCE, RYBOHR_FORCE gint allok = 0; if (g_ascii_strncasecmp("N", buff, 1) == 0) *siesta_unit_type = N_FORCE; else if (g_ascii_strncasecmp("eV/Ang", buff, 6) == 0) *siesta_unit_type = EVANG_FORCE; else if (g_ascii_strncasecmp("Ry/Bohr", buff, 7) == 0) *siesta_unit_type = RYBOHR_FORCE; else { //no match - set default *siesta_unit_type = N_FORCE; allok = 1; } return allok; } gint set_mominert_units(gchar * buff, enum siesta_mominert_measurement_type * siesta_unit_type) { //KGM_MOMINERT, RYFS_MOMINERT gint allok = 0; if (g_ascii_strncasecmp("Kg*m**2", buff, 7) == 0) *siesta_unit_type = KGM_MOMINERT; else if (g_ascii_strncasecmp("Ry*fs**2", buff, 8) == 0) *siesta_unit_type = RYFS_MOMINERT; else { //no match - set default *siesta_unit_type = KGM_MOMINERT; allok = 1; } return allok; } gint set_true_false(gchar * buff, gboolean * siesta_boolean) { gint allok = 0; if (g_ascii_strncasecmp("true", buff, 4) == 0) *siesta_boolean = TRUE; else if (g_ascii_strncasecmp(".true", buff, 5) == 0) *siesta_boolean = TRUE; else if (g_ascii_strncasecmp(".false", buff, 6) == 0) *siesta_boolean = FALSE; else if (g_ascii_strncasecmp("false", buff, 5) == 0) *siesta_boolean = FALSE; else { //no match - set default *siesta_boolean = FALSE; allok = 1; } return allok; } gint print_file_energy_units(FILE * fp, enum siesta_energy_measurement_type units) { gint allok = 0; switch (units) { case J_ENRG: fprintf(fp, " J\n"); break; case ERG_ENRG: fprintf(fp, " erg\n"); break; case EV_ENRG: fprintf(fp, " eV\n"); break; case MEV_ENRG: fprintf(fp, " meV\n"); break; case RY_ENRG: fprintf(fp, " Ry\n"); break; case MRY_ENRG: fprintf(fp, " mRy\n"); break; case HARTREE_ENRG: fprintf(fp, " Hartree\n"); break; case KO_ENRG: fprintf(fp, " K\n"); break; case KCALMOLE_ENRG: fprintf(fp, " kcal/mol\n"); break; case MHARTREE_ENRG: fprintf(fp, " mHartree\n"); break; case KJMOL_ENRG: fprintf(fp, " kJ/mol\n"); break; case HZ_ENRG: fprintf(fp, " Hz\n"); break; case THZ_ENRG: fprintf(fp, " THz\n"); break; case CM_ENRG: fprintf(fp, " cm-1\n"); break; case CMCM_ENRG: fprintf(fp, " cm**-1\n"); break; default: //no units given? //print new line fprintf(fp," \n"); allok = 1; break; } return allok; } gint print_file_mass_units(FILE * fp, enum siesta_mass_measurement_type units) { gint allok = 0; switch (units) { case KG_MASS: fprintf(fp, " Kg\n"); break; case G_MASS: fprintf(fp, " g\n"); break; case AMU_MASS: fprintf(fp, " amu\n"); break; default: //no units given? print new line fprintf(fp," \n"); allok = 1; break; } return allok; } gint print_file_length_units(FILE * fp, enum siesta_length_measurement_type units) { gint allok = 0; switch (units) //M_LEN, CM_LEN, NM_LEN, ANG_LEN, BOHR_LEN { case M_LEN: fprintf(fp, " m\n"); break; case CM_LEN: fprintf(fp, " cm\n"); break; case NM_LEN: fprintf(fp, " nm\n"); break; case ANG_LEN: fprintf(fp, " Ang\n"); break; case BOHR_LEN: fprintf(fp, " Bohr\n"); break; default: //no units given? print new line fprintf(fp," \n"); allok = 1; break; } return allok; } gint print_file_time_units(FILE * fp, enum siesta_time_measurement_type units) { gint allok = 0; switch (units) //S_TIME, FS_TIME, PS_TIME, NS_TIME { case S_TIME: fprintf(fp, " s\n"); break; case FS_TIME: fprintf(fp, " fs\n"); break; case PS_TIME: fprintf(fp, " ps\n"); break; case NS_TIME: fprintf(fp, " ns\n"); break; default: //no units given? print new line or force units? fprintf(fp," \n"); allok = 1; break; } return allok; } //Pressure gint print_file_pressure_units(FILE * fp, enum siesta_pressure_measurement_type units) { gint allok = 0; switch (units) { case P_PRES: fprintf(fp, " Pa\n"); break; case MPA_PRES: fprintf(fp, " MPa\n"); break; case GPA_PRES: fprintf(fp, " GPa\n"); break; case ATM_PRES: fprintf(fp, " atm\n"); break; case BAR_PRES: fprintf(fp, " bar\n"); break; case MBAR_PRES: fprintf(fp, " Mbar\n"); break; case RYBOHR3_PRES: fprintf(fp, " Ry/Bohr**3\n"); break; case EVANG3_PRES: fprintf(fp, " eV/Ang**3\n"); break; default: //no units given? print new line or force units? fprintf(fp," \n"); allok = 1; break; } return allok; } gint print_file_temperature_units(FILE * fp, enum siesta_temperature_measurement_type units) { gint allok = 0; switch (units) { case C_TEMP: fprintf(fp, " C\n"); break; case K_TEMP: fprintf(fp, " K\n"); break; default: //no units given? print new line or force units? fprintf(fp," \n"); allok = 1; break; } return allok; } gint print_file_force_units(FILE * fp, enum siesta_force_measurement_type units) { gint allok = 0; switch (units) //N_FORCE, EVANG_FORCE, RYBOHR_FORCE { case N_FORCE: fprintf(fp, " N\n"); break; case EVANG_FORCE: fprintf(fp, " eV/Ang\n"); break; case RYBOHR_FORCE: fprintf(fp, " Ry/Bohr\n"); break; default: //no units given? print new line or force units? fprintf(fp," \n"); allok = 1; break; } return allok; } gint print_file_mominert_units(FILE * fp, enum siesta_mominert_measurement_type units) { gint allok = 0; switch (units) //KGM_MOMINERT, RYFS_MOMINERT { case KGM_MOMINERT: fprintf(fp, " Kg*m**2\n"); break; case RYFS_MOMINERT: fprintf(fp, " Ry*fs**2\n"); break; default: //no units given? print new line or force units? fprintf(fp," \n"); allok = 1; break; } return allok; }