update fix restrain to use {fmt} and be 64-bit atom id compatible

This commit is contained in:
Axel Kohlmeyer
2020-06-03 23:43:54 -04:00
parent 9297211b06
commit aab4bc8ac2
2 changed files with 58 additions and 94 deletions

View File

@ -29,6 +29,7 @@
#include "math_const.h" #include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "fmt/format.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
@ -82,8 +83,8 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[iarg],"bond") == 0) { if (strcmp(arg[iarg],"bond") == 0) {
if (iarg+6 > narg) error->all(FLERR,"Illegal fix restrain command"); if (iarg+6 > narg) error->all(FLERR,"Illegal fix restrain command");
rstyle[nrestrain] = BOND; rstyle[nrestrain] = BOND;
ids[nrestrain][0] = force->inumeric(FLERR,arg[iarg+1]); ids[nrestrain][0] = force->tnumeric(FLERR,arg[iarg+1]);
ids[nrestrain][1] = force->inumeric(FLERR,arg[iarg+2]); ids[nrestrain][1] = force->tnumeric(FLERR,arg[iarg+2]);
kstart[nrestrain] = force->numeric(FLERR,arg[iarg+3]); kstart[nrestrain] = force->numeric(FLERR,arg[iarg+3]);
kstop[nrestrain] = force->numeric(FLERR,arg[iarg+4]); kstop[nrestrain] = force->numeric(FLERR,arg[iarg+4]);
deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]); deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]);
@ -97,8 +98,8 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"lbound") == 0) { } else if (strcmp(arg[iarg],"lbound") == 0) {
if (iarg+6 > narg) error->all(FLERR,"Illegal fix restrain command"); if (iarg+6 > narg) error->all(FLERR,"Illegal fix restrain command");
rstyle[nrestrain] = LBOUND; rstyle[nrestrain] = LBOUND;
ids[nrestrain][0] = force->inumeric(FLERR,arg[iarg+1]); ids[nrestrain][0] = force->tnumeric(FLERR,arg[iarg+1]);
ids[nrestrain][1] = force->inumeric(FLERR,arg[iarg+2]); ids[nrestrain][1] = force->tnumeric(FLERR,arg[iarg+2]);
kstart[nrestrain] = force->numeric(FLERR,arg[iarg+3]); kstart[nrestrain] = force->numeric(FLERR,arg[iarg+3]);
kstop[nrestrain] = force->numeric(FLERR,arg[iarg+4]); kstop[nrestrain] = force->numeric(FLERR,arg[iarg+4]);
deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]); deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]);
@ -112,9 +113,9 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"angle") == 0) { } else if (strcmp(arg[iarg],"angle") == 0) {
if (iarg+7 > narg) error->all(FLERR,"Illegal fix restrain command"); if (iarg+7 > narg) error->all(FLERR,"Illegal fix restrain command");
rstyle[nrestrain] = ANGLE; rstyle[nrestrain] = ANGLE;
ids[nrestrain][0] = force->inumeric(FLERR,arg[iarg+1]); ids[nrestrain][0] = force->tnumeric(FLERR,arg[iarg+1]);
ids[nrestrain][1] = force->inumeric(FLERR,arg[iarg+2]); ids[nrestrain][1] = force->tnumeric(FLERR,arg[iarg+2]);
ids[nrestrain][2] = force->inumeric(FLERR,arg[iarg+3]); ids[nrestrain][2] = force->tnumeric(FLERR,arg[iarg+3]);
kstart[nrestrain] = force->numeric(FLERR,arg[iarg+4]); kstart[nrestrain] = force->numeric(FLERR,arg[iarg+4]);
kstop[nrestrain] = force->numeric(FLERR,arg[iarg+5]); kstop[nrestrain] = force->numeric(FLERR,arg[iarg+5]);
target[nrestrain] = force->numeric(FLERR,arg[iarg+6]); target[nrestrain] = force->numeric(FLERR,arg[iarg+6]);
@ -124,10 +125,10 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
if (iarg+8 > narg) error->all(FLERR,"Illegal fix restrain command"); if (iarg+8 > narg) error->all(FLERR,"Illegal fix restrain command");
rstyle[nrestrain] = DIHEDRAL; rstyle[nrestrain] = DIHEDRAL;
mult[nrestrain] = 1; mult[nrestrain] = 1;
ids[nrestrain][0] = force->inumeric(FLERR,arg[iarg+1]); ids[nrestrain][0] = force->tnumeric(FLERR,arg[iarg+1]);
ids[nrestrain][1] = force->inumeric(FLERR,arg[iarg+2]); ids[nrestrain][1] = force->tnumeric(FLERR,arg[iarg+2]);
ids[nrestrain][2] = force->inumeric(FLERR,arg[iarg+3]); ids[nrestrain][2] = force->tnumeric(FLERR,arg[iarg+3]);
ids[nrestrain][3] = force->inumeric(FLERR,arg[iarg+4]); ids[nrestrain][3] = force->tnumeric(FLERR,arg[iarg+4]);
kstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]); kstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]);
kstop[nrestrain] = force->numeric(FLERR,arg[iarg+6]); kstop[nrestrain] = force->numeric(FLERR,arg[iarg+6]);
target[nrestrain] = force->numeric(FLERR,arg[iarg+7]); target[nrestrain] = force->numeric(FLERR,arg[iarg+7]);
@ -271,24 +272,16 @@ void FixRestrain::restrain_bond(int m)
if (newton_bond) { if (newton_bond) {
if (i2 == -1 || i2 >= nlocal) return; if (i2 == -1 || i2 >= nlocal) return;
if (i1 == -1) { if (i1 == -1)
char str[128]; error->one(FLERR,fmt::format("Restrain atoms {} {} missing on "
sprintf(str, "proc {} at step {}", ids[m][0],ids[m][1],
"Restrain atoms %d %d missing on proc %d at step " BIGINT_FORMAT, comm->me,update->ntimestep));
ids[m][0],ids[m][1],
comm->me,update->ntimestep);
error->one(FLERR,str);
}
} else { } else {
if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal)) return; if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal)) return;
if (i1 == -1 || i2 == -1) { if (i1 == -1 || i2 == -1)
char str[128]; error->one(FLERR,fmt::format("Restrain atoms {} {} missing on "
sprintf(str, "proc {} at step {}", ids[m][0],ids[m][1],
"Restrain atoms %d %d missing on proc %d at step " BIGINT_FORMAT, comm->me,update->ntimestep));
ids[m][0],ids[m][1],
comm->me,update->ntimestep);
error->one(FLERR,str);
}
} }
delx = x[i1][0] - x[i2][0]; delx = x[i1][0] - x[i2][0];
@ -352,24 +345,16 @@ void FixRestrain::restrain_lbound(int m)
if (newton_bond) { if (newton_bond) {
if (i2 == -1 || i2 >= nlocal) return; if (i2 == -1 || i2 >= nlocal) return;
if (i1 == -1) { if (i1 == -1)
char str[128]; error->one(FLERR,fmt::format("Restrain atoms {} {} missing on "
sprintf(str, "proc {} at step {}",ids[m][0],ids[m][1],
"Restrain atoms %d %d missing on proc %d at step " BIGINT_FORMAT, comm->me,update->ntimestep));
ids[m][0],ids[m][1],
comm->me,update->ntimestep);
error->one(FLERR,str);
}
} else { } else {
if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal)) return; if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal)) return;
if (i1 == -1 || i2 == -1) { if (i1 == -1 || i2 == -1)
char str[128]; error->one(FLERR,fmt::format("Restrain atoms {} {} missing on "
sprintf(str, "proc {} at step {}",ids[m][0],ids[m][1],
"Restrain atoms %d %d missing on proc %d at step " BIGINT_FORMAT, comm->me,update->ntimestep));
ids[m][0],ids[m][1],
comm->me,update->ntimestep);
error->one(FLERR,str);
}
} }
delx = x[i1][0] - x[i2][0]; delx = x[i1][0] - x[i2][0];
@ -442,27 +427,17 @@ void FixRestrain::restrain_angle(int m)
if (newton_bond) { if (newton_bond) {
if (i2 == -1 || i2 >= nlocal) return; if (i2 == -1 || i2 >= nlocal) return;
if (i1 == -1 || i3 == -1) { if (i1 == -1 || i3 == -1)
char str[128]; error->one(FLERR,fmt::format("Restrain atoms {} {} {} missing on "
sprintf(str, "proc {} at step {}",ids[m][0],ids[m][1],
"Restrain atoms %d %d %d missing on proc %d at step " ids[m][2],comm->me,update->ntimestep));
BIGINT_FORMAT,
ids[m][0],ids[m][1],ids[m][2],
comm->me,update->ntimestep);
error->one(FLERR,str);
}
} else { } else {
if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal) && if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal) &&
(i3 == -1 || i3 >= nlocal)) return; (i3 == -1 || i3 >= nlocal)) return;
if (i1 == -1 || i2 == -1 || i3 == -1) { if (i1 == -1 || i2 == -1 || i3 == -1)
char str[128]; error->one(FLERR,fmt::format("Restrain atoms {} {} {} missing on "
sprintf(str, "proc {} at step {}",ids[m][0],ids[m][1],
"Restrain atoms %d %d %d missing on proc %d at step " ids[m][2],comm->me,update->ntimestep));
BIGINT_FORMAT,
ids[m][0],ids[m][1],ids[m][2],
comm->me,update->ntimestep);
error->one(FLERR,str);
}
} }
// 1st bond // 1st bond
@ -572,27 +547,19 @@ void FixRestrain::restrain_dihedral(int m)
if (newton_bond) { if (newton_bond) {
if (i2 == -1 || i2 >= nlocal) return; if (i2 == -1 || i2 >= nlocal) return;
if (i1 == -1 || i3 == -1 || i4 == -1) { if (i1 == -1 || i3 == -1 || i4 == -1)
char str[128]; error->one(FLERR,fmt::format("Restrain atoms {} {} {} {} missing on "
sprintf(str, "proc {} at step {}",ids[m][0],ids[m][1],
"Restrain atoms %d %d %d %d missing on proc %d at step " ids[m][2],ids[m][3],comm->me,
BIGINT_FORMAT, update->ntimestep));
ids[m][0],ids[m][1],ids[m][2],ids[m][3],
comm->me,update->ntimestep);
error->one(FLERR,str);
}
} else { } else {
if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal) && if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal) &&
(i3 == -1 || i3 >= nlocal) && (i4 == -1 || i3 >= nlocal)) return; (i3 == -1 || i3 >= nlocal) && (i4 == -1 || i3 >= nlocal)) return;
if (i1 == -1 || i2 == -1 || i3 == -1 || i4 == -1) { if (i1 == -1 || i2 == -1 || i3 == -1 || i4 == -1)
char str[128]; error->one(FLERR,fmt::format("Restrain atoms {} {} {} {} missing on "
sprintf(str, "proc {} at step {}",ids[m][0],ids[m][1],
"Restrain atoms %d %d %d %d missing on proc %d at step " ids[m][2],ids[m][3],comm->me,
BIGINT_FORMAT, update->ntimestep));
ids[m][0],ids[m][1],ids[m][2],ids[m][3],
comm->me,update->ntimestep);
error->one(FLERR,str);
}
} }
// 1st bond // 1st bond
@ -648,21 +615,18 @@ void FixRestrain::restrain_dihedral(int m)
int me; int me;
MPI_Comm_rank(world,&me); MPI_Comm_rank(world,&me);
if (screen) { if (screen) {
char str[128]; std::string str = fmt::format("Restrain problem: {} {} {} {} {} {}",
sprintf(str,"Restrain problem: %d " BIGINT_FORMAT " " me,update->ntimestep,atom->tag[i1],
TAGINT_FORMAT " " TAGINT_FORMAT " " atom->tag[i2],atom->tag[i3],atom->tag[i4]);
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str); error->warning(FLERR,str);
fprintf(screen," 1st atom: %d %g %g %g\n", fmt::print(screen," 1st atom: {} {} {} {}\n"
me,x[i1][0],x[i1][1],x[i1][2]); " 2nd atom: {} {} {} {}\n"
fprintf(screen," 2nd atom: %d %g %g %g\n", " 3rd atom: {} {} {} {}\n"
me,x[i2][0],x[i2][1],x[i2][2]); " 4th atom: {} {} {} {}\n",
fprintf(screen," 3rd atom: %d %g %g %g\n", comm->me,x[i1][0],x[i1][1],x[i1][2],
me,x[i3][0],x[i3][1],x[i3][2]); comm->me,x[i2][0],x[i2][1],x[i2][2],
fprintf(screen," 4th atom: %d %g %g %g\n", comm->me,x[i3][0],x[i3][1],x[i3][2],
me,x[i4][0],x[i4][1],x[i4][2]); comm->me,x[i4][0],x[i4][1],x[i4][2]);
} }
} }

View File

@ -43,7 +43,7 @@ class FixRestrain : public Fix {
int nrestrain,maxrestrain; int nrestrain,maxrestrain;
int *rstyle; int *rstyle;
int *mult; int *mult;
int **ids; tagint **ids;
double *kstart,*kstop,*deqstart,*deqstop,*target; double *kstart,*kstop,*deqstart,*deqstop,*target;
double *cos_target,*sin_target; double *cos_target,*sin_target;
double energy,energy_all; double energy,energy_all;