diff --git a/src/SMTBQ/pair_smtbq.cpp b/src/SMTBQ/pair_smtbq.cpp index be8af7f9bd..f5649b9214 100644 --- a/src/SMTBQ/pair_smtbq.cpp +++ b/src/SMTBQ/pair_smtbq.cpp @@ -51,6 +51,7 @@ #include "memory.h" #include "neigh_list.h" #include "neighbor.h" +#include "text_file_reader.h" #include "update.h" #include @@ -73,6 +74,8 @@ using namespace MathSpecial; #define PGDELTA 1 #define MAXNEIGH 24 +static constexpr char SMTBQ_SEPARATORS[] = "' \t\n\r"; + /* ---------------------------------------------------------------------- */ PairSMTBQ::PairSMTBQ(LAMMPS *lmp) : Pair(lmp) @@ -130,6 +133,12 @@ PairSMTBQ::PairSMTBQ(LAMMPS *lmp) : Pair(lmp) Nevery = 0.0; Neverypot = 0.0; + QEqMode = nullptr; + Bavard = nullptr; + QInitMode = nullptr; + writepot = nullptr; + writeenerg = nullptr; + fct = nullptr; maxpage = 0; @@ -148,16 +157,16 @@ PairSMTBQ::~PairSMTBQ() { int i; if (elements) { - for (i = 0; i < atom->ntypes ; i++ ) free( params[i].nom); - for (i = 1; i <= maxintparam ; i++ ) free( intparams[i].typepot); - for (i = 1; i <= maxintparam ; i++ ) free( intparams[i].mode); + for (i = 0; i < atom->ntypes ; i++ ) delete[] params[i].nom; + for (i = 1; i <= maxintparam ; i++ ) delete[] intparams[i].typepot; + for (i = 1; i <= maxintparam ; i++ ) delete[] intparams[i].mode; } - free(QEqMode); - free(QInitMode); - free(writepot); - free(writeenerg); - free(Bavard); + delete[] QEqMode; + delete[] Bavard; + delete[] QInitMode; + delete[] writepot; + delete[] writeenerg; memory->sfree(params); memory->sfree(intparams); @@ -307,7 +316,17 @@ double PairSMTBQ::init_one(int i, int j) void PairSMTBQ::read_file(char *file) { int num_atom_types,i,k,m,test,j,verbose; - char **words; + + if (params) { + for (i = 0; i < atom->ntypes ; i++ ) delete[] params[i].nom; + for (i = 1; i <= maxintparam ; i++ ) delete[] intparams[i].typepot; + for (i = 1; i <= maxintparam ; i++ ) delete[] intparams[i].mode; + delete[] QEqMode; + delete[] Bavard; + delete[] QInitMode; + delete[] writepot; + delete[] writeenerg; + } memory->sfree(params); params = nullptr; @@ -325,444 +344,462 @@ void PairSMTBQ::read_file(char *file) coordOxSurf = 0.0; ROxBB = 0.0; ROxSurf = 0.0; + QEqMode = nullptr; + Bavard = nullptr; + QInitMode = nullptr; + writepot = nullptr; + writeenerg = nullptr; - // open file on all processors - FILE *fp; - fp = utils::open_potential(file,lmp,nullptr); - if (fp == nullptr) { - char str[128]; - snprintf(str,128,"Cannot open SMTBQ potential file %s",file); - error->one(FLERR,str); + // open file on all processors + auto potfile = utils::get_potential_file_path(file); + if (!potfile.empty()) { + auto date = utils::get_potential_date(potfile, "SMTBQ"); + auto units = utils::get_potential_units(potfile, "SMTBQ"); + if (!date.empty() && (comm->me == 0)) + utils::logmesg(lmp, "Reader SMTBQ potential file {} with DATE: {}\n",potfile,date); + if (!units.empty() && (units != update->unit_style) && (comm->me == 0)) + error->one(FLERR, "Potential file {} requires {} units but {} units are in use", + potfile, units, update->unit_style); + } else { + error->all(FLERR, "Cannot open SMTBQ potential file {}", file); } + TextFileReader reader(potfile, "SMTBQ"); - // read each line out of file, skipping blank lines or leading '#' - // store line of params if all 3 element tags are in element list + try { - char *ptr; + // Nombre d'atome different dans la structure + // =============================================== + const char *line = reader.next_line(); + auto values = ValueTokenizer(line, SMTBQ_SEPARATORS); + auto label = values.next_string(); + num_atom_types = values.next_int(); + if ((comm->me == 0) && verbose) utils::logmesg(lmp," {} {}\n",label,num_atom_types); - ptr = (char*) malloc(sizeof(char)*MAXLINE); - words = (char**) malloc(sizeof(char*)*MAXTOKENS); - for (i=0; i < MAXTOKENS; i++) - words[i] = (char*) malloc(sizeof(char)*MAXTOKENS); - - - /* strip comment, skip line if blank */ - - if (verbose) printf ("\n"); - fgets(ptr,MAXLINE,fp); - while (strchr(ptr,'#')) { - if (verbose) printf ("%s",ptr); - fgets(ptr,MAXLINE,fp); - } - - - // Nombre d'atome different dans la structure - // =============================================== - Tokenize( ptr, &words ); - num_atom_types = atoi(words[1]); - if (verbose) printf (" %s %d\n", words[0], num_atom_types); - - memory->create(intype,num_atom_types,num_atom_types,"pair:intype"); - - m = 0; - for (i = 0; i < num_atom_types; i++) { - for (j = 0; j < num_atom_types; j++) { - if (j < i) { intype[i][j] = intype[j][i];} - else { intype[i][j] = 0; - m = m + 1; } - if (verbose) printf ("i %d, j %d, intype %d - nb pair %d\n",i,j,intype[i][j],m); - } - } - - // load up parameter settings and error check their values - - nparams = maxparam = num_atom_types; - params = (Param *) memory->create(params,maxparam*sizeof(Param), - "pair:params"); - maxintparam = m; - intparams = (Intparam *) memory->create(intparams,(maxintparam+1)*sizeof(Intparam), - "pair:intparams"); - - for (i=0; i < num_atom_types; i++) - params[i].nom = (char*) malloc(sizeof(char)*3); - - for (i=1; i <= maxintparam; i++) - intparams[i].typepot = (char*) malloc(sizeof(char)*15); - - for (i=1; i <= maxintparam; i++) - intparams[i].mode = (char*) malloc(sizeof(char)*6); - - QEqMode = (char*) malloc(sizeof(char)*19); - Bavard = (char*) malloc(sizeof(char)*6); - QInitMode = (char*) malloc(sizeof(char)*19); - writepot = (char*) malloc(sizeof(char)*6); - writeenerg = (char*) malloc(sizeof(char)*6); - - - // Little loop for ion's parameters - // ================================================ - for (i=0; icreate(intype,num_atom_types,num_atom_types,"pair:intype"); + m = 0; + for (i = 0; i < num_atom_types; i++) { + for (j = 0; j < num_atom_types; j++) { + if (j < i) { intype[i][j] = intype[j][i]; } + else { intype[i][j] = 0; m = m + 1; } + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, "i {}, j {}, intype {} - nb pair {}\n",i,j,intype[i][j],m); + } } + // load up parameter settings and error check their values + nparams = maxparam = num_atom_types; + maxintparam = m; + params = (Param *) memory->create(params,maxparam*sizeof(Param),"pair:params"); + intparams = (Intparam *) memory->create(intparams,(maxintparam+1)*sizeof(Intparam), + "pair:intparams"); - // Line 4bis - Coordinance et rayon pour Ox - if (strcmp(params[i].nom,"O")==0) { + for (i=0; i < num_atom_types; i++) + params[i].nom = nullptr; - fgets( ptr, MAXLINE, fp); - Tokenize ( ptr, &words ); + for (i=1; i <= maxintparam; i++) + intparams[i].typepot = nullptr; - coordOxBB= atof(words[1]) ; - coordOxBulk= atof(words[2]) ; - coordOxSurf= atof(words[3]) ; - ROxBB = atof(words[4]) ; - params[i].R = atof(words[5]) ; - ROxSurf = atof(words[6]) ; - if (verbose) printf(" %s %f %f %f %f %f %f\n",words[0],coordOxBB,coordOxBulk,coordOxSurf,ROxBB,params[i].R,ROxSurf); - } + for (i=1; i <= maxintparam; i++) + intparams[i].mode = nullptr; - // Ligne 5 - Nombre d'etats partages - fgets( ptr, MAXLINE, fp); - Tokenize ( ptr, &words ); - params[i].n0 = atof(words[1]); - if (verbose) printf(" %s %f\n",words[0],params[i].n0); + QEqMode = nullptr; + Bavard = nullptr; + QInitMode = nullptr; + writepot = nullptr; + writeenerg = nullptr; - // Parametres de Slater - params[i].dzeta = (2.0*params[i].ne + 1.0)/(4.0*params[i].R); - if (verbose) printf (" Parametre dzeta (Slater) : %f\n",params[i].dzeta); - } // Fin elements i + // Little loop for ion's parameters + // ================================================ + for (i=0; ime == 0) && verbose) utils::logmesg(lmp, label+"\n"); - /* ===================================================================== - reading the interaction's parameters - ===================================================================== */ + // Line 2 - Al - m = 0; maxintsm = 0; // - for (k=0 ; k<=maxintparam ; k++) {intparams[k].intsm = 0;} - // --------------------------------- - for (k = 0; k < maxintparam; k++) { + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + params[i].nom = utils::strdup(values.next_string()); + params[i].sto = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n", label,params[i].nom,params[i].sto); + + //Line 3 - Charges + + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + params[i].qform = values.next_double(); + params[i].masse = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n", label,params[i].qform,params[i].masse); + + // Line 4 - Parametres QEq + + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + params[i].ne = values.next_double(); + params[i].chi = values.next_double(); + params[i].dj = values.next_double(); + + if (strcmp(params[i].nom, "O") !=0) { + params[i].R = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {} {} {}\n",label,params[i].ne,params[i].chi, + params[i].dj,params[i].R); + } else { + if ((comm->me == 0) && verbose) + utils::logmesg(lmp," {} {} {} {}\n",label,params[i].ne,params[i].chi,params[i].dj); + } + + + // Line 4bis - Coordinance et rayon pour Ox + if (strcmp(params[i].nom,"O") == 0) { + + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + coordOxBB = values.next_double(); + coordOxBulk = values.next_double(); + coordOxSurf = values.next_double(); + ROxBB = values.next_double(); + params[i].R = values.next_double(); + ROxSurf = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {} {} {} {} {}\n", + label,coordOxBB,coordOxBulk,coordOxSurf,ROxBB,params[i].R,ROxSurf); + } + + // Ligne 5 - Nombre d'etats partages + + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + params[i].n0 = values.next_double(); + if ((comm->me == 0) && verbose) utils::logmesg(lmp, " {} {}\n",label,params[i].n0); + + // Parametres de Slater + params[i].dzeta = (2.0*params[i].ne + 1.0)/(4.0*params[i].R); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " Parametre dzeta (Slater) : {}\n",params[i].dzeta); + + } // Fin elements i + + /* ===================================================================== + reading the interaction's parameters + ===================================================================== */ + + m = 0; maxintsm = 0; // + for (k=0 ; k<=maxintparam ; k++) {intparams[k].intsm = 0;} // --------------------------------- - m += 1; + for (k = 0; k < maxintparam; k++) { + // --------------------------------- + m += 1; - // Ligne 5 - parametre des potentials - fgets(ptr,MAXLINE,fp); if (verbose) printf ("%s",ptr); + // Ligne 5 - parametre des potentials + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + if ((comm->me == 0) && verbose) utils::logmesg(lmp, label+"\n"); - // Lecture des protagonistes - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); + // Lecture des protagonistes + test = 0; + auto words = Tokenizer(reader.next_line(), SMTBQ_SEPARATORS).as_vector(); + if (words.size() < 3) test = 1; - test = 0; - for (i = 0; i %d = %s\n",words[1],i,params[i].nom); - for (j = 0; j %d = %s\n",words[2],j,params[j].nom); - - if (test == 1) { - if (verbose) printf ("========== fin des interaction ==========\n"); - break ; } - - - intype[i][j] = m; - intype[j][i] = intype[i][j]; - strcpy( intparams[m].typepot , words[3] ); - intparams[m].intsm = 0; - if (verbose) printf (" itype %d jtype %d - intype %d\n",i,j,intype[i][j]); - - if (strcmp(intparams[m].typepot,"second_moment") !=0 && - strcmp(intparams[m].typepot,"buck") != 0 && - strcmp(intparams[m].typepot,"buckPlusAttr") != 0) { - error->all(FLERR,"the potential other than second_moment or buckingham have not treated here\n");} - - - // On detemrine le type d'interaction - // ----------------------------------- - if (strcmp(intparams[m].typepot,"second_moment") == 0) { - maxintsm += 1; - strcpy( intparams[m].mode , words[4] ); - intparams[m].intsm = maxintsm; - - if (strcmp(intparams[m].mode,"oxide") != 0 && - strcmp(intparams[m].mode,"metal") != 0) { - error->all(FLERR,"needs mode to second moment interaction : oxide or metal"); } - - // if (strcmp(intparams[m].mode,"oxide") == 0) - // intparams[m].ncov = min((params[i].sto)*(params[i].n0),(params[j].sto)*(params[j].n0)); - - if (verbose) printf(" %s %s %s %s %s \n",words[0],words[1],words[2], - intparams[m].typepot,intparams[m].mode); - - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); - - intparams[m].a = atof(words[1]) ; - intparams[m].p = atof(words[2]) ; - intparams[m].ksi = atof(words[3]) ; - intparams[m].q = atof(words[4]) ; - if (verbose) printf (" %s %f %f %f %f\n",words[0], - intparams[m].a,intparams[m].p,intparams[m].ksi,intparams[m].q); - - // Ligne 6 - rayon de coupure potential SM - - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); - - intparams[m].dc1 = atof(words[1]) ; - intparams[m].dc2 = atof(words[2]) ; - intparams[m].r0 = atof(words[3]) ; - - - if (strcmp(intparams[m].mode,"metal") == 0) { - if (verbose) printf (" %s %f %f %f\n",words[0], - intparams[m].dc1,intparams[m].dc2,intparams[m].r0); - } else { - if (verbose) printf (" %s %f %f %f\n",words[0], - intparams[m].dc1,intparams[m].dc2,intparams[m].r0); + if (test == 1) { + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, "========== fin des interaction ==========\n"); + break ; } + intype[i][j] = m; + intype[j][i] = intype[i][j]; + intparams[m].typepot = utils::strdup(words[3]); + intparams[m].intsm = 0; + if ((comm->me == 0) && verbose) + utils::logmesg(lmp," itype {} jtype {} - intype {}\n",i,j,intype[i][j]); - } else if (strcmp(intparams[m].typepot,"buck") == 0) { + if (strcmp(intparams[m].typepot,"second_moment") != 0 && + strcmp(intparams[m].typepot,"buck") != 0 && + strcmp(intparams[m].typepot,"buckPlusAttr") != 0) { + error->all(FLERR,"Potential must be 'second_moment' or 'buckingham'\n"); + } - if (verbose) printf(" %s %s %s %s\n",words[0],words[1],words[2], - intparams[m].typepot); + // On detemrine le type d'interaction + // ----------------------------------- + if (strcmp(intparams[m].typepot,"second_moment") == 0) { + maxintsm += 1; + intparams[m].mode = utils::strdup(words[4]); + intparams[m].intsm = maxintsm; - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); + if (strcmp(intparams[m].mode,"oxide") != 0 && + strcmp(intparams[m].mode,"metal") != 0) { + error->all(FLERR,"Mode of second moment interaction must be 'oxide' or 'metal'"); + } - intparams[m].abuck = atof(words[1]) ; intparams[m].rhobuck = atof(words[2]) ; - if (verbose) printf (" %s %f %f\n",words[0],intparams[m].abuck,intparams[m].rhobuck); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {} {} {}\n",words[0],words[1],words[2], + intparams[m].typepot,intparams[m].mode); + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + intparams[m].a = values.next_double(); + intparams[m].p = values.next_double(); + intparams[m].ksi = values.next_double(); + intparams[m].q = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {} {} {}\n",label,intparams[m].a, + intparams[m].p,intparams[m].ksi,intparams[m].q); + + // Ligne 6 - rayon de coupure potential SM + + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + intparams[m].dc1 = values.next_double(); + intparams[m].dc2 = values.next_double(); + intparams[m].r0 = values.next_double(); + + + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {} {}\n",label, intparams[m].dc1, + intparams[m].dc2,intparams[m].r0); + + } else if (strcmp(intparams[m].typepot,"buck") == 0) { + + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {} {}\n",words[0],words[1],words[2],intparams[m].typepot); + + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + intparams[m].abuck = values.next_double(); + intparams[m].rhobuck = values.next_double(); + + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n",label,intparams[m].abuck,intparams[m].rhobuck); + + } else if (strcmp(intparams[m].typepot,"buckPlusAttr") == 0) { + + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {} {}\n",words[0],words[1],words[2],intparams[m].typepot); + + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + intparams[m].abuck = values.next_double(); + intparams[m].rhobuck = values.next_double(); + + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n",label,intparams[m].abuck,intparams[m].rhobuck); + + line = reader.next_line(); + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + intparams[m].aOO = values.next_double(); + intparams[m].bOO = values.next_double(); + intparams[m].r1OO = values.next_double(); + intparams[m].r2OO = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {} {} {}\n",label,intparams[m].aOO, + intparams[m].bOO,intparams[m].r1OO,intparams[m].r2OO); + + } + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " intsm {} \n",intparams[m].intsm); + + } // for maxintparam + + /* ==================================================================== + tables Parameters + ==================================================================== */ + // Ligne 9 - rayon de coupure Electrostatique + if (test == 0) { + line = reader.next_line(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, line); + line = reader.next_line(); } - else if (strcmp(intparams[m].typepot,"buckPlusAttr") == 0) { + values = ValueTokenizer(line, SMTBQ_SEPARATORS); + label = values.next_string(); + cutmax = values.next_double(); + for (i=0 ; ime == 0) && verbose) + utils::logmesg(lmp, " {} {} \n",label,cutmax); - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); + // Ligne 9 - parametre pour les tableaux - intparams[m].abuck = atof(words[1]) ; intparams[m].rhobuck = atof(words[2]) ; - if (verbose) printf (" %s %f %f\n",words[0],intparams[m].abuck,intparams[m].rhobuck); + values = reader.next_values(0,SMTBQ_SEPARATORS); + label = values.next_string(); + rmin = values.next_double(); + dr = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n",label, rmin, dr); + kmax = int(cutmax*cutmax/(2.0*dr*rmin)); + ds = cutmax*cutmax/static_cast(kmax) ; + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " kmax {} et ds {}\n",kmax,ds); - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); + /* ======================================================== */ + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, reader.next_line()); + else reader.skip_line(); - intparams[m].aOO = atof(words[1]) ; intparams[m].bOO = atof(words[2]) ; - intparams[m].r1OO = atof(words[3]) ;intparams[m].r2OO = atof(words[4]) ; - if (verbose) printf (" %s %f %f %f %f \n",words[0],intparams[m].aOO, - intparams[m].bOO,intparams[m].r1OO,intparams[m].r2OO); + values = reader.next_values(0, SMTBQ_SEPARATORS); + label = values.next_string(); + Qstep = values.next_bigint(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {}\n",label, Qstep); + values = reader.next_values(0, SMTBQ_SEPARATORS); + label = values.next_string(); + loopmax = values.next_int(); + precision = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n",label, loopmax, precision); + + /* Param de coordination ============================================= */ + + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, reader.next_values(0, SMTBQ_SEPARATORS).next_string()+"\n"); + else reader.skip_line(); + + values = reader.next_values(0, SMTBQ_SEPARATORS); + label = values.next_string(); + r1Coord = values.next_double(); + r2Coord = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n",label,r1Coord,r2Coord); + + /* Mode for QInit============================================= */ + + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, reader.next_values(0, SMTBQ_SEPARATORS).next_string()+"\n"); + else reader.skip_line(); + + values = reader.next_values(0, SMTBQ_SEPARATORS); + label = values.next_string(); + QInitMode = utils::strdup(values.next_string()); + if (strcmp(QInitMode,"true") == 0) QOxInit= values.next_double(); + else QOxInit = 0.0; + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n",label,QInitMode,QOxInit); + + /* Mode for QEq============================================= */ + + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, reader.next_values(0, SMTBQ_SEPARATORS).next_string()+"\n"); + else reader.skip_line(); + + values = reader.next_values(0, SMTBQ_SEPARATORS); + label = values.next_string(); + QEqMode = utils::strdup(values.next_string()); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {}\n",label,QEqMode); + + if (strcmp(QEqMode,"BulkFromSlab") == 0) { + zlim1QEq = values.next_double(); + zlim2QEq = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n",label,zlim1QEq,zlim2QEq); + + } else if (strcmp(QEqMode,"Surface") == 0) { + + zlim1QEq = values.next_double(); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {}\n",label,zlim1QEq); + + } else if (strcmp(QEqMode,"QEqAll") != 0 && + strcmp(QEqMode,"QEqAllParallel") != 0 && + strcmp(QEqMode,"Surface") != 0) { + error->all(FLERR,"The QEq Mode is not known. QEq mode should be :\n" + " Possible QEq modes | parameters\n" + " QEqAll | no parameters\n" + " QEqAllParallel | no parameters\n" + " Surface | zlim (QEq only for z>zlim)\n" + " BulkFromSlab | zlim1 zlim2 (QEq only for zlim1me == 0) && verbose) + utils::logmesg(lmp, reader.next_values(0, SMTBQ_SEPARATORS).next_string()+"\n"); + else reader.skip_line(); - /* ==================================================================== - tables Parameters - ==================================================================== */ + values = reader.next_values(0, SMTBQ_SEPARATORS); + label = values.next_string(); + Bavard = utils::strdup(values.next_string()); + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {}\n",label,Bavard); - // Ligne 9 - rayon de coupure Electrostatique - if (test == 0) { - fgets(ptr,MAXLINE,fp); - if (verbose) printf ("%s\n",ptr); + // --------------------------------------- + // Writing the energy component. - fgets( ptr, MAXLINE, fp); - } - Tokenize( ptr, &words ); + values = reader.next_values(0, SMTBQ_SEPARATORS); + label = values.next_string(); + writeenerg = utils::strdup(values.next_string()); + if (strcmp(writeenerg,"true") == 0) Nevery = values.next_double(); + else Nevery = 0.0; + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n",label,writeenerg,Nevery); - for (i=0 ; ime == 0) && verbose) + utils::logmesg(lmp, " {} {} {}\n",label,writepot,Neverypot); - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); + // === Rayon de coupure premier voisins : 1,2*r0 + for (i=0 ; i(kmax) ; - if (verbose) printf (" kmax %d et ds %f\n",kmax,ds); - - /* ======================================================== */ - fgets( ptr, MAXLINE, fp); - if (verbose) printf ("%s",ptr); - - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); - Qstep = atoi(words[1]); - if (verbose) printf (" %s " BIGINT_FORMAT "\n",words[0],Qstep); - - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); - loopmax = atoi(words[1]); - precision = atof(words[2]); - if (verbose) printf (" %s %d %f\n",words[0],loopmax,precision); - - /* Param de coordination ============================================= */ - - fgets( ptr, MAXLINE, fp); - if (verbose) printf ("%s",ptr); - - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); - r1Coord = atof(words[1]); - r2Coord = atof(words[2]); - if (verbose) printf (" %s %f %f\n",words[0],r1Coord,r2Coord); - - /* Mode for QInit============================================= */ - fgets( ptr, MAXLINE, fp); - if (verbose) printf ("%s",ptr); - - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); - strcpy( QInitMode , words[1] ); - if (strcmp(QInitMode,"true") == 0) QOxInit= atof(words[2]); - else QOxInit = 0.0; - if (verbose) printf (" %s %s %f\n",words[0],QInitMode,QOxInit); - - - /* Mode for QEq============================================= */ - - fgets( ptr, MAXLINE, fp); - if (verbose) printf ("%s",ptr); - - fgets( ptr, MAXLINE, fp); - Tokenize( ptr, &words ); - strcpy( QEqMode , words[1] ); - if (verbose) printf (" %s %s\n",words[0],QEqMode); - - fgets( ptr, MAXLINE, fp); - - if (strcmp(QEqMode,"BulkFromSlab") == 0) { - Tokenize( ptr, &words ); - zlim1QEq = atof(words[1]); - zlim2QEq = atof(words[2]); - if (verbose) printf (" %s %f %f\n",words[0],zlim1QEq,zlim2QEq); - - } else if (strcmp(QEqMode,"Surface") == 0) { - Tokenize( ptr, &words ); - zlim1QEq = atof(words[1]); - if (verbose) printf (" %s %f \n",words[0],zlim1QEq); - - } else if (strcmp(QEqMode,"QEqAll") != 0 && - strcmp(QEqMode,"QEqAllParallel") != 0 && - strcmp(QEqMode,"Surface") != 0) { - error->all(FLERR,"The QEq Mode is not known. QEq mode should be :\n" - " Possible QEq modes | parameters\n" - " QEqAll | no parameters\n" - " QEqAllParallel | no parameters\n" - " Surface | zlim (QEq only for z>zlim)\n" - " BulkFromSlab | zlim1 zlim2 (QEq only for zlim1 %f Ang\n", - intparams[m].typepot,intparams[m].neig_cut); + intparams[m].neig_cut = 1.2*intparams[m].r0; + if (strcmp(intparams[m].typepot,"second_moment") == 0 ) + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " Rc 1er voisin, typepot {} -> {} Ang\n", + intparams[m].typepot,intparams[m].neig_cut); + } } + + //A adapter au STO + for (i=1,ncov=params[0].sto*params[0].n0; i < nparams; ++i) + ncov = min(ncov,(params[1].sto)*(params[1].n0)); + + if ((comm->me == 0) && verbose) + utils::logmesg(lmp, " Parametre ncov = {}\n" + " ********************************************* \n",ncov); + } catch (std::exception &e) { + error->all(FLERR, "Error parsing SMTBQ parameter file: {}", e.what()); } - - //A adapter au STO - for (i=1,ncov=params[0].sto*params[0].n0; i < nparams; ++i) - ncov = min(ncov,(params[1].sto)*(params[1].n0)); - - if (verbose) printf (" Parametre ncov = %f\n",ncov); - if (verbose) printf (" ********************************************* \n"); } /* ---------------------------------------------------------------------- @@ -1113,7 +1150,7 @@ void PairSMTBQ::compute(int eflag, int vflag) if (me == 0&& strcmp(Bavard,"false") != 0) { printf ("A la fin de Compute\n"); - printf ("Nemton_pair : %d, evflag %d, tail_flag %d,vflag_fdotr %d\n", + printf ("Newton_pair : %d, evflag %d, tail_flag %d,vflag_fdotr %d\n", newton_pair,evflag,force->pair->tail_flag,vflag_fdotr); printf ("neighbor->includegroup %d\n",neighbor->includegroup); @@ -3473,28 +3510,6 @@ void PairSMTBQ::add_pages(int howmany) /* ---------------------------------------------------------------------- */ - -int PairSMTBQ::Tokenize( char* s, char*** tok ) -{ - char test[MAXLINE]; - const char *sep = "' "; - char *mot; - int count=0; - mot = nullptr; - - - strncpy( test, s, MAXLINE-1 ); - - for (mot = strtok(test, sep); mot; mot = strtok(nullptr, sep)) { - strncpy( (*tok)[count], mot, MAXLINE ); - count++; - } - - return count; -} - - - void PairSMTBQ::CheckEnergyVSForce() { double drL,iq,jq,rsq,evdwlCoul,fpairCoul,eflag=0,ErepR,frepR,fpair,evdwl; diff --git a/src/SMTBQ/pair_smtbq.h b/src/SMTBQ/pair_smtbq.h index 636cff8266..2809405697 100644 --- a/src/SMTBQ/pair_smtbq.h +++ b/src/SMTBQ/pair_smtbq.h @@ -149,8 +149,6 @@ class PairSMTBQ : public Pair { void forward_int(int *); void reverse_int(int *); - int Tokenize(char *, char ***); - template const T &min(const T &a, const T &b) { return !(b < a) ? a : b; // or: return !comp(b,a)?a:b; for the comp version