diff --git a/examples/SPIN/gneb/iron/final.iron_spin b/examples/SPIN/gneb/iron/final.iron_spin index aa1cbae770..97aca57199 100644 --- a/examples/SPIN/gneb/iron/final.iron_spin +++ b/examples/SPIN/gneb/iron/final.iron_spin @@ -32,37 +32,3 @@ 31 2.2000000000000002e+00 7.1662499999999998e+00 1.0032750000000000e+01 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 32 2.2000000000000002e+00 1.0032750000000000e+01 1.0032750000000000e+01 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/SPIN/gneb/iron/initial.iron_spin b/examples/SPIN/gneb/iron/initial.iron_spin index 1d28afe400..30f4fa39f0 100644 --- a/examples/SPIN/gneb/iron/initial.iron_spin +++ b/examples/SPIN/gneb/iron/initial.iron_spin @@ -5,7 +5,7 @@ LAMMPS data file via write_data, version 4 Jan 2019, timestep = 0 0.0000000000000000e+00 1.1465999999999999e+01 xlo xhi 0.0000000000000000e+00 1.1465999999999999e+01 ylo yhi -0.0000000000000000e+00 2.8664999999999998e+00 zlo zhi +-1.0000000000000000e+00 10.0664999999999998e+00 zlo zhi Masses diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 8910959355..bb15e99fd2 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -135,24 +135,24 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : else procnext = -1; uworld = universe->uworld; - int *iroots = new int[nreplica]; - MPI_Group uworldgroup,rootgroup; if (NEBLongRange) { - for (int i=0; iroot_proc[i]; + int *iroots = new int[nreplica]; + MPI_Group uworldgroup,rootgroup; + + for (int i=0; iroot_proc[i]; MPI_Comm_group(uworld, &uworldgroup); MPI_Group_incl(uworldgroup, nreplica, iroots, &rootgroup); MPI_Comm_create(uworld, rootgroup, &rootworld); + if (rootgroup != MPI_GROUP_NULL) MPI_Group_free(&rootgroup); + if (uworldgroup != MPI_GROUP_NULL) MPI_Group_free(&uworldgroup); + delete [] iroots; } - delete [] iroots; // create a new compute pe style // id = fix-ID + pe, compute group = all - std::string cmd = id + std::string("_pe"); - id_pe = new char[cmd.size()+1]; - strcpy(id_pe,cmd.c_str()); - modify->add_compute(cmd + " all pe"); + id_pe = utils::strdup(std::string(id)+"_pe"); + modify->add_compute(std::string(id_pe)+" all pe"); // initialize local storage diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index e730aac351..3136bdae7c 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -27,6 +27,7 @@ #include "output.h" #include "thermo.h" #include "timer.h" +#include "tokenizer.h" #include "universe.h" #include "update.h" @@ -173,12 +174,11 @@ void NEB::run() else color = 1; MPI_Comm_split(uworld,color,0,&roots); - int ineb; - for (ineb = 0; ineb < modify->nfix; ineb++) - if (strcmp(modify->fix[ineb]->style,"neb") == 0) break; - if (ineb == modify->nfix) error->all(FLERR,"NEB requires use of fix neb"); + auto fixes = modify->get_fix_by_style("^neb$"); + if (fixes.size() != 1) + error->all(FLERR,"NEB requires use of exactly one fix neb instance"); - fneb = (FixNEB *) modify->fix[ineb]; + fneb = (FixNEB *) fixes[0]; if (verbose) numall =7; else numall = 4; memory->create(all,nreplica,numall,"neb:all"); @@ -376,11 +376,11 @@ void NEB::run() void NEB::readfile(char *file, int flag) { - int i,j,m,nchunk,eofflag,nlines; + int i,nchunk,eofflag,nlines; tagint tag; char *eof,*start,*next,*buf; char line[MAXLINE]; - double xx,yy,zz,delx,dely,delz; + double delx,dely,delz; if (me_universe == 0 && universe->uscreen) fprintf(universe->uscreen,"Reading NEB coordinate file(s) ...\n"); @@ -424,10 +424,7 @@ void NEB::readfile(char *file, int flag) } char *buffer = new char[CHUNK*MAXLINE]; - char **values = new char*[ATTRIBUTE_PERLINE]; - double fraction = ireplica/(nreplica-1.0); - double **x = atom->x; int nlocal = atom->nlocal; @@ -435,9 +432,7 @@ void NEB::readfile(char *file, int flag) // two versions of read_lines_from_file() for world vs universe bcast // count # of atom coords changed so can check for invalid atom IDs in file - int ncount = 0; - - int nread = 0; + int ncount = 0, nread = 0; while (nread < nlines) { nchunk = MIN(nlines-nread,CHUNK); if (flag == 0) @@ -461,49 +456,47 @@ void NEB::readfile(char *file, int flag) for (i = 0; i < nchunk; i++) { next = strchr(buf,'\n'); + *next = '\0'; - values[0] = strtok(buf," \t\n\r\f"); - for (j = 1; j < nwords; j++) - values[j] = strtok(nullptr," \t\n\r\f"); + try { + ValueTokenizer values(buf," \t\n\r\f"); - // adjust atom coord based on replica fraction - // for flag = 0, interpolate for intermediate and final replicas - // for flag = 1, replace existing coord with new coord - // ignore image flags of replica x - // displacement from first replica is via minimum image convention - // if x of some replica is across periodic boundary: - // new x may be outside box - // will be remapped back into box when simulation starts - // its image flags will then be adjusted + // adjust atom coord based on replica fraction + // for flag = 0, interpolate for intermediate and final replicas + // for flag = 1, replace existing coord with new coord + // ignore image flags of replica x + // displacement from first replica is via minimum image convention + // if x of some replica is across periodic boundary: + // new x may be outside box + // will be remapped back into box when simulation starts + // its image flags will then be adjusted - tag = ATOTAGINT(values[0]); - m = atom->map(tag); - if (m >= 0 && m < nlocal) { - ncount++; - xx = atof(values[1]); - yy = atof(values[2]); - zz = atof(values[3]); + tag = values.next_tagint(); + int m = atom->map(tag); + if (m >= 0 && m < nlocal) { + ncount++; - delx = xx - x[m][0]; - dely = yy - x[m][1]; - delz = zz - x[m][2]; + delx = values.next_double() - x[m][0]; + dely = values.next_double() - x[m][1]; + delz = values.next_double() - x[m][2]; - domain->minimum_image(delx,dely,delz); + domain->minimum_image(delx,dely,delz); - if (flag == 0) { - x[m][0] += fraction*delx; - x[m][1] += fraction*dely; - x[m][2] += fraction*delz; - } else { - x[m][0] += delx; - x[m][1] += dely; - x[m][2] += delz; + if (flag == 0) { + x[m][0] += fraction*delx; + x[m][1] += fraction*dely; + x[m][2] += fraction*delz; + } else { + x[m][0] += delx; + x[m][1] += dely; + x[m][2] += delz; + } } + } catch (std::exception &e) { + error->universe_one(FLERR,"Incorrectly formatted NEB file: " + std::string(e.what())); } - buf = next + 1; } - nread += nchunk; } @@ -522,9 +515,7 @@ void NEB::readfile(char *file, int flag) } // clean up - - delete [] buffer; - delete [] values; + delete[] buffer; if (flag == 0) { if (me_universe == 0) { diff --git a/src/SPIN/fix_neb_spin.cpp b/src/SPIN/fix_neb_spin.cpp index de50cac764..0a96227c46 100644 --- a/src/SPIN/fix_neb_spin.cpp +++ b/src/SPIN/fix_neb_spin.cpp @@ -56,7 +56,7 @@ FixNEBSpin::FixNEBSpin(LAMMPS *lmp, int narg, char **arg) : counts(nullptr), displacements(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal fix neb_spin command"); + if (narg < 4) error->all(FLERR,"Illegal fix neb/spin command"); kspring = utils::numeric(FLERR,arg[3],false,lmp); if (kspring <= 0.0) error->all(FLERR,"Illegal fix neb command"); @@ -77,16 +77,16 @@ FixNEBSpin::FixNEBSpin(LAMMPS *lmp, int narg, char **arg) : int iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"parallel") == 0) { - error->all(FLERR,"Illegal fix neb command"); + error->all(FLERR,"Illegal fix neb/spin command"); iarg += 2; } else if (strcmp(arg[iarg],"perp") == 0) { - error->all(FLERR,"Illegal fix neb command"); + error->all(FLERR,"Illegal fix neb/spin command"); iarg += 2; } else if (strcmp (arg[iarg],"end") == 0) { iarg += 3; } else if (strcmp (arg[iarg],"lattice") == 0) { iarg += 2; - } else error->all(FLERR,"Illegal fix neb command"); + } else error->all(FLERR,"Illegal fix neb/spin command"); } // nreplica = number of partitions @@ -107,16 +107,18 @@ FixNEBSpin::FixNEBSpin(LAMMPS *lmp, int narg, char **arg) : else procnext = -1; uworld = universe->uworld; - int *iroots = new int[nreplica]; - MPI_Group uworldgroup,rootgroup; if (NEBLongRange) { - for (int i=0; iroot_proc[i]; + int *iroots = new int[nreplica]; + MPI_Group uworldgroup,rootgroup; + + for (int i=0; iroot_proc[i]; MPI_Comm_group(uworld, &uworldgroup); MPI_Group_incl(uworldgroup, nreplica, iroots, &rootgroup); MPI_Comm_create(uworld, rootgroup, &rootworld); + if (rootgroup != MPI_GROUP_NULL) MPI_Group_free(&rootgroup); + if (uworldgroup != MPI_GROUP_NULL) MPI_Group_free(&uworldgroup); + delete [] iroots; } - delete [] iroots; // create a new compute pe style // id = fix-ID + pe, compute group = all diff --git a/src/SPIN/neb_spin.cpp b/src/SPIN/neb_spin.cpp index 0ed3ad28c3..8c375e632b 100644 --- a/src/SPIN/neb_spin.cpp +++ b/src/SPIN/neb_spin.cpp @@ -38,6 +38,7 @@ #include "output.h" #include "thermo.h" #include "timer.h" +#include "tokenizer.h" #include "universe.h" #include "update.h" @@ -147,6 +148,7 @@ void NEBSpin::command(int narg, char **arg) verbose=false; if (strcmp(arg[narg-1],"verbose") == 0) verbose=true; + // run the NEB calculation run(); } @@ -166,12 +168,11 @@ void NEBSpin::run() // search for neb_spin fix, allocate it - int ineb; - for (ineb = 0; ineb < modify->nfix; ineb++) - if (strcmp(modify->fix[ineb]->style,"neb/spin") == 0) break; - if (ineb == modify->nfix) error->all(FLERR,"NEBSpin requires use of fix neb/spin"); + auto fixes = modify->get_fix_by_style("^neb/spin"); + if (fixes.size() != 1) + error->all(FLERR,"NEBSpin requires use of exactly one fix neb/spin instance"); - fneb = (FixNEBSpin *) modify->fix[ineb]; + fneb = (FixNEBSpin *) fixes[0]; if (verbose) numall =7; else numall = 4; memory->create(all,nreplica,numall,"neb:all"); @@ -370,12 +371,11 @@ void NEBSpin::run() void NEBSpin::readfile(char *file, int flag) { - int i,j,m,nchunk,eofflag,nlines; + int i,nchunk,eofflag,nlines; tagint tag; char *eof,*start,*next,*buf; char line[MAXLINE]; - double xx,yy,zz; - double musp,spx,spy,spz; + double musp,xx,yy,zz,spx,spy,spz; if (me_universe == 0 && universe->uscreen) fprintf(universe->uscreen,"Reading NEBSpin coordinate file(s) ...\n"); @@ -393,10 +393,12 @@ void NEBSpin::readfile(char *file, int flag) start = &line[strspn(line," \t\n\v\f\r")]; if (*start != '\0' && *start != '#') break; } - sscanf(line,"%d",&nlines); + int rv = sscanf(line,"%d",&nlines); + if (rv != 1) nlines = -1; } MPI_Bcast(&nlines,1,MPI_INT,0,uworld); - + if (nlines < 0) + error->universe_all(FLERR,"Incorrectly formatted NEB file"); } else { if (me == 0) { if (ireplica) { @@ -407,17 +409,17 @@ void NEBSpin::readfile(char *file, int flag) start = &line[strspn(line," \t\n\v\f\r")]; if (*start != '\0' && *start != '#') break; } - sscanf(line,"%d",&nlines); + int rv = sscanf(line,"%d",&nlines); + if (rv != 1) nlines = -1; } else nlines = 0; } MPI_Bcast(&nlines,1,MPI_INT,0,world); + if (nlines < 0) + error->all(FLERR,"Incorrectly formatted NEB file"); } char *buffer = new char[CHUNK*MAXLINE]; - char **values = new char*[ATTRIBUTE_PERLINE]; - double fraction = ireplica/(nreplica-1.0); - double **x = atom->x; double **sp = atom->sp; double spinit[3],spfinal[3]; @@ -427,11 +429,7 @@ void NEBSpin::readfile(char *file, int flag) // two versions of read_lines_from_file() for world vs universe bcast // count # of atom coords changed so can check for invalid atom IDs in file - int ncount = 0; - - int temp_flag,rot_flag; - temp_flag = rot_flag = 0; - int nread = 0; + int ncount=0, nread=0, temp_flag=0, rot_flag=0; while (nread < nlines) { nchunk = MIN(nlines-nread,CHUNK); if (flag == 0) @@ -444,7 +442,7 @@ void NEBSpin::readfile(char *file, int flag) buf = buffer; next = strchr(buf,'\n'); *next = '\0'; - int nwords = utils::trim_and_count_words(buf); + int nwords = utils::count_words(utils::trim_comment(buf)); *next = '\n'; if (nwords != ATTRIBUTE_PERLINE) @@ -455,75 +453,77 @@ void NEBSpin::readfile(char *file, int flag) for (i = 0; i < nchunk; i++) { next = strchr(buf,'\n'); + *next = '\0'; - values[0] = strtok(buf," \t\n\r\f"); - for (j = 1; j < nwords; j++) - values[j] = strtok(nullptr," \t\n\r\f"); + try { + ValueTokenizer values(buf," \t\n\r\f"); - // adjust spin coord based on replica fraction - // for flag = 0, interpolate for intermediate and final replicas - // for flag = 1, replace existing coord with new coord - // ignore image flags of final x - // for interpolation: - // new x is displacement from old x via minimum image convention - // if final x is across periodic boundary: - // new x may be outside box - // will be remapped back into box when simulation starts - // its image flags will then be adjusted + // adjust spin coord based on replica fraction + // for flag = 0, interpolate for intermediate and final replicas + // for flag = 1, replace existing coord with new coord + // ignore image flags of final x + // for interpolation: + // new x is displacement from old x via minimum image convention + // if final x is across periodic boundary: + // new x may be outside box + // will be remapped back into box when simulation starts + // its image flags will then be adjusted - tag = ATOTAGINT(values[0]); - m = atom->map(tag); - if (m >= 0 && m < nlocal) { - ncount++; - musp = atof(values[1]); - xx = atof(values[2]); - yy = atof(values[3]); - zz = atof(values[4]); - spx = atof(values[5]); - spy = atof(values[6]); - spz = atof(values[7]); + tag = values.next_tagint(); + int m = atom->map(tag); + if (m >= 0 && m < nlocal) { + ncount++; - if (flag == 0) { + musp = values.next_double(); + xx = values.next_double(); + yy = values.next_double(); + zz = values.next_double(); + spx = values.next_double(); + spy = values.next_double(); + spz = values.next_double(); - spinit[0] = sp[m][0]; - spinit[1] = sp[m][1]; - spinit[2] = sp[m][2]; - spfinal[0] = spx; - spfinal[1] = spy; - spfinal[2] = spz; + if (flag == 0) { - // interpolate intermediate spin states + spinit[0] = sp[m][0]; + spinit[1] = sp[m][1]; + spinit[2] = sp[m][2]; + spfinal[0] = spx; + spfinal[1] = spy; + spfinal[2] = spz; - sp[m][3] = musp; - if (fraction == 0.0) { - sp[m][0] = spinit[0]; - sp[m][1] = spinit[1]; - sp[m][2] = spinit[2]; - } else if (fraction == 1.0) { - sp[m][0] = spfinal[0]; - sp[m][1] = spfinal[1]; - sp[m][2] = spfinal[2]; + // interpolate intermediate spin states + + sp[m][3] = musp; + if (fraction == 0.0) { + sp[m][0] = spinit[0]; + sp[m][1] = spinit[1]; + sp[m][2] = spinit[2]; + } else if (fraction == 1.0) { + sp[m][0] = spfinal[0]; + sp[m][1] = spfinal[1]; + sp[m][2] = spfinal[2]; + } else { + temp_flag = initial_rotation(spinit,spfinal,fraction); + rot_flag = MAX(temp_flag,rot_flag); + sp[m][0] = spfinal[0]; + sp[m][1] = spfinal[1]; + sp[m][2] = spfinal[2]; + } } else { - temp_flag = initial_rotation(spinit,spfinal,fraction); - rot_flag = MAX(temp_flag,rot_flag); - sp[m][0] = spfinal[0]; - sp[m][1] = spfinal[1]; - sp[m][2] = spfinal[2]; + sp[m][3] = musp; + x[m][0] = xx; + x[m][1] = yy; + x[m][2] = zz; + sp[m][0] = spx; + sp[m][1] = spy; + sp[m][2] = spz; } - } else { - sp[m][3] = musp; - x[m][0] = xx; - x[m][1] = yy; - x[m][2] = zz; - sp[m][0] = spx; - sp[m][1] = spy; - sp[m][2] = spz; } + } catch (std::exception &e) { + error->universe_one(FLERR,"Incorrectly formatted NEB file: " + std::string(e.what())); } - buf = next + 1; } - nread += nchunk; } @@ -549,9 +549,7 @@ void NEBSpin::readfile(char *file, int flag) } // clean up - delete[] buffer; - delete[] values; if (flag == 0) { if (me_universe == 0) {