diff --git a/doc/src/Commands_all.rst b/doc/src/Commands_all.rst index 92c0ba107a..458238cca6 100644 --- a/doc/src/Commands_all.rst +++ b/doc/src/Commands_all.rst @@ -94,6 +94,7 @@ An alphabetic list of all general LAMMPS commands. * :doc:`package ` * :doc:`pair_coeff ` * :doc:`pair_modify ` + * :doc:`pair_style ` * :doc:`pair_write ` * :doc:`partition ` * :doc:`prd ` diff --git a/src/KIM/kim_init.cpp b/src/KIM/kim_init.cpp index 5328769d22..9f16a2bcdb 100644 --- a/src/KIM/kim_init.cpp +++ b/src/KIM/kim_init.cpp @@ -57,6 +57,8 @@ ------------------------------------------------------------------------- */ #include "kim_init.h" +#include "fix_store_kim.h" +#include "kim_units.h" #include #include #include @@ -71,8 +73,7 @@ #include "input.h" #include "variable.h" #include "citeme.h" -#include "fix_store_kim.h" -#include "kim_units.h" +#include "utils.h" extern "C" { #include "KIM_SimulatorHeaders.h" diff --git a/src/USER-PLUMED/fix_plumed.cpp b/src/USER-PLUMED/fix_plumed.cpp index e1f0ea0bfe..3a52bbaaa1 100644 --- a/src/USER-PLUMED/fix_plumed.cpp +++ b/src/USER-PLUMED/fix_plumed.cpp @@ -82,6 +82,7 @@ FixPlumed::FixPlumed(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Incompatible API version for PLUMED in fix plumed. " "Only Plumed 2.4.x, 2.5.x, and 2.6.x are tested and supported."); +#if !defined(MPI_STUBS) // If the -partition option is activated then enable // inter-partition communication @@ -108,7 +109,6 @@ FixPlumed::FixPlumed(LAMMPS *lmp, int narg, char **arg) : // whereas if partitions are not defined then world is equal to // MPI_COMM_WORLD. -#if !defined(MPI_STUBS) // plumed does not know about LAMMPS using the MPI STUBS library and will // fail if this is called under these circumstances p->cmd("setMPIComm",&world); diff --git a/src/USER-QTB/fix_qbmsst.cpp b/src/USER-QTB/fix_qbmsst.cpp index 5ac611b13e..d58c917308 100644 --- a/src/USER-QTB/fix_qbmsst.cpp +++ b/src/USER-QTB/fix_qbmsst.cpp @@ -21,6 +21,7 @@ #include #include #include +#include #include "atom.h" #include "force.h" #include "update.h" @@ -34,6 +35,7 @@ #include "kspace.h" #include "math_const.h" #include "utils.h" +#include "fmt/format.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -98,58 +100,71 @@ FixQBMSST::FixQBMSST(LAMMPS *lmp, int narg, char **arg) : while (iarg < narg) { if (strcmp(arg[iarg],"q") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - qmass = atof(arg[iarg+1]); if (qmass < 0.0) error->all(FLERR,"Fix qbmsst qmass must be >= 0.0"); + qmass = force->numeric(FLERR,arg[iarg+1]); + if (qmass < 0.0) error->all(FLERR,"Fix qbmsst qmass must be >= 0.0"); iarg += 2; } else if (strcmp(arg[iarg],"mu") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - mu = atof(arg[iarg+1]); if (mu < 0.0) error->all(FLERR,"Fix qbmsst mu must be >= 0.0"); + mu = force->numeric(FLERR,arg[iarg+1]); + if (mu < 0.0) error->all(FLERR,"Fix qbmsst mu must be >= 0.0"); iarg += 2; } else if (strcmp(arg[iarg],"p0") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - p0 = atof(arg[iarg+1]); if (p0 < 0.0) error->all(FLERR,"Fix qbmsst p0 must be >= 0.0"); + p0 = force->numeric(FLERR,arg[iarg+1]); + if (p0 < 0.0) error->all(FLERR,"Fix qbmsst p0 must be >= 0.0"); p0_set = 1; iarg += 2; } else if (strcmp(arg[iarg],"v0") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - v0 = atof(arg[iarg+1]); if (v0 < 0.0) error->all(FLERR,"Fix qbmsst v0 must be >= 0.0"); + v0 = force->numeric(FLERR,arg[iarg+1]); + if (v0 < 0.0) error->all(FLERR,"Fix qbmsst v0 must be >= 0.0"); v0_set = 1; iarg += 2; } else if (strcmp(arg[iarg],"e0") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - e0 = atof(arg[iarg+1]); + e0 = force->numeric(FLERR,arg[iarg+1]); e0_set = 1; iarg += 2; } else if (strcmp(arg[iarg],"tscale") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - tscale = atof(arg[iarg+1]); if (tscale < 0.0 || tscale > 1.0) error->all(FLERR,"Fix qbmsst tscale must satisfy 0 <= tscale < 1"); + tscale = force->numeric(FLERR,arg[iarg+1]); + if (tscale < 0.0 || tscale > 1.0) error->all(FLERR,"Fix qbmsst tscale must satisfy 0 <= tscale < 1"); iarg += 2; } else if (strcmp(arg[iarg],"damp") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - t_period = atof(arg[iarg+1]); if (t_period <= 0.0) error->all(FLERR,"Fix qbmsst damp must be > 0.0"); fric_coef = 1/t_period; + t_period = force->numeric(FLERR,arg[iarg+1]); + if (t_period <= 0.0) error->all(FLERR,"Fix qbmsst damp must be > 0.0"); + fric_coef = 1/t_period; iarg += 2; } else if (strcmp(arg[iarg],"seed") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - seed = atof(arg[iarg+1]); if (seed <= 0) error->all(FLERR,"Fix qbmsst seed must be a positive integer"); + seed = force->inumeric(FLERR,arg[iarg+1]); + if (seed <= 0) error->all(FLERR,"Fix qbmsst seed must be a positive integer"); iarg += 2; } else if (strcmp(arg[iarg],"f_max") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - f_max = atof(arg[iarg+1]); if (f_max <= 0) error->all(FLERR,"Fix qbmsst f_max must be > 0.0"); + f_max = force->numeric(FLERR,arg[iarg+1]); + if (f_max <= 0) error->all(FLERR,"Fix qbmsst f_max must be > 0.0"); iarg += 2; } else if (strcmp(arg[iarg],"N_f") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - N_f = atof(arg[iarg+1]); if (N_f <= 0) error->all(FLERR,"Fix qbmsst N_f must be a positive integer"); + N_f = force->inumeric(FLERR,arg[iarg+1]); + if (N_f <= 0) error->all(FLERR,"Fix qbmsst N_f must be a positive integer"); iarg += 2; } else if (strcmp(arg[iarg],"eta") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - eta = atof(arg[iarg+1]); if (eta <= 0) error->all(FLERR,"Fix qbmsst eta must be >= 0.0"); + eta = force->numeric(FLERR,arg[iarg+1]); + if (eta <= 0) error->all(FLERR,"Fix qbmsst eta must be >= 0.0"); iarg += 2; } else if (strcmp(arg[iarg],"beta") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - beta = atof(arg[iarg+1]); if (beta <= 0) error->all(FLERR,"Fix qbmsst beta must be a positive integer"); + beta = force->inumeric(FLERR,arg[iarg+1]); + if (beta <= 0) error->all(FLERR,"Fix qbmsst beta must be a positive integer"); iarg += 2; } else if (strcmp(arg[iarg],"T_init") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qbmsst command"); - t_init = atof(arg[iarg+1]); if (t_init <= 0) error->all(FLERR,"Fix qbmsst T_init must be >= 0.0"); + t_init = force->numeric(FLERR,arg[iarg+1]); + if (t_init <= 0) error->all(FLERR,"Fix qbmsst T_init must be >= 0.0"); iarg += 2; } else error->all(FLERR,"Illegal fix qbmsst command"); } @@ -157,48 +172,32 @@ FixQBMSST::FixQBMSST(LAMMPS *lmp, int narg, char **arg) : // check for periodicity in controlled dimensions if (domain->nonperiodic) error->all(FLERR,"Fix qbmsst requires a periodic box"); + maxexchange = 6*N_f+3; + // comments if (comm->me == 0) { - if (screen) { - fprintf(screen,"QBMSST parameters:\n"); - if (direction == 0) fprintf(screen," Shock in x direction\n"); - else if (direction == 1) fprintf(screen," Shock in y direction\n"); - else if (direction == 2) fprintf(screen," Shock in z direction\n"); + std::string msg = "QBMSST parameters:\n"; - fprintf(screen," Cell mass-like parameter qmass (units of mass^2/length^4) = %12.5e\n", qmass); - fprintf(screen," Shock velocity = %12.5e\n", velocity); - fprintf(screen," Artificial viscosity (units of mass/length/time) = %12.5e\n", mu); + if (direction == 0) msg += " Shock in x direction\n"; + else if (direction == 1) msg += " Shock in y direction\n"; + else if (direction == 2) msg += " Shock in z direction\n"; - if (p0_set) - fprintf(screen," Initial pressure specified to be %12.5e\n", p0); - else fprintf(screen," Initial pressure calculated on first step\n"); - if (v0_set) - fprintf(screen," Initial volume specified to be %12.5e\n", v0); - else fprintf(screen," Initial volume calculated on first step\n"); - if (e0_set) - fprintf(screen," Initial energy specified to be %12.5e\n", e0); - else fprintf(screen," Initial energy calculated on first step\n"); - } - if (logfile) { - fprintf(logfile,"QBMSST parameters:\n"); - if (direction == 0) fprintf(logfile," Shock in x direction\n"); - else if (direction == 1) fprintf(logfile," Shock in y direction\n"); - else if (direction == 2) fprintf(logfile," Shock in z direction\n"); + msg += fmt::format(" Cell mass-like parameter qmass " + "(units of mass^2/length^4) = {:12.5e}\n", qmass); + msg += fmt::format(" Shock velocity = {:12.5e}\n", velocity); + msg += fmt::format(" Artificial viscosity (units of " + "mass/length/time) = {:12.5e}\n", mu); - fprintf(logfile," Cell mass-like parameter qmass (units of mass^2/length^4) = %12.5e\n", qmass); - fprintf(logfile," Shock velocity = %12.5e\n", velocity); - fprintf(logfile," Artificial viscosity (units of mass/length/time) = %12.5e\n", mu); - - if (p0_set) - fprintf(logfile," Initial pressure specified to be %12.5e\n", p0); - else fprintf(logfile," Initial pressure calculated on first step\n"); - if (v0_set) - fprintf(logfile," Initial volume specified to be %12.5e\n", v0); - else fprintf(logfile," Initial volume calculated on first step\n"); - if (e0_set) - fprintf(logfile," Initial energy specified to be %12.5e\n", e0); - else fprintf(logfile," Initial energy calculated on first step\n"); - } + if (p0_set) + msg += fmt::format(" Initial pressure specified to be {:12.5e}\n", p0); + else msg += " Initial pressure calculated on first step\n"; + if (v0_set) + msg += fmt::format(" Initial volume specified to be {:12.5e}\n", v0); + else msg += " Initial volume calculated on first step\n"; + if (e0_set) + msg += fmt::format(" Initial energy specified to be {:12.5e}\n", e0); + else msg += " Initial energy calculated on first step\n"; + utils::logmesg(lmp,msg); } // create a new compute temp style @@ -363,14 +362,17 @@ void FixQBMSST::init() //set up the h time step for updating the random force \delta{}h=\frac{\pi}{\Omega_{max}} if (int(1.0/(2*f_max*dtv)) == 0) { - if (comm->me == 0) printf ("Warning: Either f_max is too high or the time step is too big, setting f_max to be 1/timestep!\n"); + if (comm->me == 0) error->warning(FLERR,"Either f_max is too high or the time step " + "is too big, setting f_max to be 1/timestep!\n"); h_timestep=dtv; alpha=1; } else { alpha=int(1.0/(2*f_max*dtv)); h_timestep=alpha*dtv; } - if (comm->me == 0) printf ("The effective maximum frequncy is now %f inverse time unit with alpha value as %d!\n", 0.5/h_timestep, alpha); + if (comm->me == 0 && screen) + fmt::print(screen,"The effective maximum frequency is now {} inverse time unit " + "with alpha value as {}!\n", 0.5/h_timestep, alpha); //gfactor is the random force \sqrt{\frac{2\gamma{}m_{i}}{\alpha*\delta{}t}}, \sqrt{12} makes the random array variance equal to unit. for (int i = 1; i <= atom->ntypes; i++) { @@ -432,20 +434,16 @@ void FixQBMSST::setup(int /*vflag*/) if ( v0_set == 0 ) { v0 = compute_vol(); v0_set = 1; - if (comm->me == 0) { - if ( screen ) fprintf(screen,"Fix QBMSST v0 = %12.5e\n", v0); - if ( logfile ) fprintf(logfile,"Fix QBMSST v0 = %12.5e\n", v0); - } + if (comm->me == 0) + utils::logmesg(lmp,fmt::format("Fix QBMSST v0 = {:12.5e}\n", v0)); } if ( p0_set == 0 ) { p0 = p_current[direction]; p0_set = 1; - if ( comm->me == 0 ) { - if ( screen ) fprintf(screen,"Fix QBMSST p0 = %12.5e\n", p0); - if ( logfile ) fprintf(logfile,"Fix QBMSST p0 = %12.5e\n", p0); - } + if ( comm->me == 0 ) + utils::logmesg(lmp,fmt::format("Fix QBMSST p0 = {:12.5e}\n", p0)); } if ( e0_set == 0 ) { @@ -453,11 +451,8 @@ void FixQBMSST::setup(int /*vflag*/) e0_set = 1; old_eavg = e0; - if ( comm->me == 0 ) { - if ( screen ) fprintf(screen,"Fix QBMSST e0 = to be %12.5e\n",e0); - if ( logfile ) fprintf(logfile,"Fix QBMSST e0 = to be %12.5e\n",e0); - } - + if ( comm->me == 0 ) + utils::logmesg(lmp,fmt::format("Fix QBMSST e0 = to be {:12.5e}\n",e0)); } temperature->compute_vector(); @@ -476,16 +471,10 @@ void FixQBMSST::setup(int /*vflag*/) omega[direction]=-1*sqrt(fac1); double fac2 = omega[direction]/v0; - if ( comm->me == 0 && tscale != 1.0) { - if ( screen ) - fprintf(screen,"Fix QBMSST initial strain rate of %12.5e established " - "by reducing temperature by factor of %12.5e\n", - fac2,tscale); - if ( logfile ) - fprintf(logfile,"Fix QBMSST initial strain rate of %12.5e established " - "by reducing temperature by factor of %12.5e\n", - fac2,tscale); - } + if ( comm->me == 0 && tscale != 1.0) + utils::logmesg(lmp,fmt::format("Fix QBMSST initial strain rate of {:12.5e} " + "established by reducing temperature by " + "factor of {:12.5e}\n",fac2,tscale)); for (int i = 0; i < atom->nlocal; i++) { if (mask[i] & groupbit) { for (int k = 0; k < 3; k++ ) { @@ -531,30 +520,31 @@ void FixQBMSST::initial_integrate(int /*vflag*/) // decide if the qtb temperature need to be updated or not if (counter_l == 0) { t_current -= dtv*fric_coef*eta*beta*(old_eavg-e0)/(3*ntotal*boltz); - old_eavg = 0;//clear old energy average - if (t_current < 0.0) t_current=0; + if (t_current > 0.0) { + old_eavg = 0;//clear old energy average - // load omega_H with calculated spectrum at a specific temperature (corrected spectrum), omega_H is the Fourier transformation of time_H - for (int k = 0; k < 2*N_f; k++) { - double f_k=(k-N_f)/(2*N_f*h_timestep); //\omega_k=\frac{2\pi}{\delta{}h}\frac{k}{2N_f} for k from -N_f to N_f-1 - if(k == N_f) { - omega_H[k]=sqrt(force->boltz * t_current); - } else { - double energy_k= force->hplanck * fabs(f_k); - omega_H[k]=sqrt( energy_k * (0.5+1.0/( exp(energy_k/(force->boltz * t_current)) - 1.0 )) ); - omega_H[k]*=alpha*sin((k-N_f)*MY_PI/(2*alpha*N_f))/sin((k-N_f)*MY_PI/(2*N_f)); - } - } - - // construct the signal filter H, filter has the unit of of sqrt(energy) \sqrt{2N_f}^{-1}H\left(t_n\right) - for (int n = 0; n < 2*N_f; n++) { - time_H[n] = 0; - double t_n=(n-N_f); + // load omega_H with calculated spectrum at a specific temperature (corrected spectrum), omega_H is the Fourier transformation of time_H for (int k = 0; k < 2*N_f; k++) { - double omega_k=(k-N_f)*MY_PI/N_f; - time_H[n] += omega_H[k]*(cos(omega_k*t_n)); + double f_k=(k-N_f)/(2*N_f*h_timestep); //\omega_k=\frac{2\pi}{\delta{}h}\frac{k}{2N_f} for k from -N_f to N_f-1 + if(k == N_f) { + omega_H[k]=sqrt(force->boltz * t_current); + } else { + double energy_k= force->hplanck * fabs(f_k); + omega_H[k]=sqrt( energy_k * (0.5+1.0/( exp(energy_k/(force->boltz * t_current)) - 1.0 )) ); + omega_H[k]*=alpha*sin((k-N_f)*MY_PI/(2*alpha*N_f))/sin((k-N_f)*MY_PI/(2*N_f)); + } + } + + // construct the signal filter H, filter has the unit of of sqrt(energy) \sqrt{2N_f}^{-1}H\left(t_n\right) + for (int n = 0; n < 2*N_f; n++) { + time_H[n] = 0; + double t_n=(n-N_f); + for (int k = 0; k < 2*N_f; k++) { + double omega_k=(k-N_f)*MY_PI/N_f; + time_H[n] += omega_H[k]*(cos(omega_k*t_n)); + } + time_H[n]/=(2.0*N_f); } - time_H[n]/=(2.0*N_f); } } @@ -1170,11 +1160,12 @@ void FixQBMSST::copy_arrays(int i, int j, int /*delflag*/) ------------------------------------------------------------------------- */ int FixQBMSST::pack_exchange(int i, double *buf) { - for (int m = 0; m < 2*N_f; m++) buf[m] = random_array_0[i][m]; - for (int m = 0; m < 2*N_f; m++) buf[m+2*N_f] = random_array_1[i][m]; - for (int m = 0; m < 2*N_f; m++) buf[m+4*N_f] = random_array_2[i][m]; - for (int m = 0; m < 3; m++) buf[m+6*N_f] = fran[i][m]; - return 6*N_f+3; + int n = 0; + for (int m = 0; m < 2*N_f; m++) buf[n++] = random_array_0[i][m]; + for (int m = 0; m < 2*N_f; m++) buf[n++] = random_array_1[i][m]; + for (int m = 0; m < 2*N_f; m++) buf[n++] = random_array_2[i][m]; + for (int m = 0; m < 3; m++) buf[n++] = fran[i][m]; + return n; } /* ---------------------------------------------------------------------- @@ -1182,9 +1173,10 @@ int FixQBMSST::pack_exchange(int i, double *buf) ------------------------------------------------------------------------- */ int FixQBMSST::unpack_exchange(int nlocal, double *buf) { - for (int m = 0; m < 2*N_f; m++) random_array_0[nlocal][m] = buf[m]; - for (int m = 0; m < 2*N_f; m++) random_array_1[nlocal][m] = buf[m+2*N_f]; - for (int m = 0; m < 2*N_f; m++) random_array_2[nlocal][m] = buf[m+4*N_f]; - for (int m = 0; m < 3; m++) fran[nlocal][m] = buf[m+6*N_f]; - return 6*N_f+3; + int n = 0; + for (int m = 0; m < 2*N_f; m++) random_array_0[nlocal][m] = buf[n++]; + for (int m = 0; m < 2*N_f; m++) random_array_1[nlocal][m] = buf[n++]; + for (int m = 0; m < 2*N_f; m++) random_array_2[nlocal][m] = buf[n++]; + for (int m = 0; m < 3; m++) fran[nlocal][m] = buf[n++]; + return n; } diff --git a/src/USER-QTB/fix_qtb.cpp b/src/USER-QTB/fix_qtb.cpp index fa15385859..646cf06ce2 100644 --- a/src/USER-QTB/fix_qtb.cpp +++ b/src/USER-QTB/fix_qtb.cpp @@ -32,6 +32,8 @@ #include "math_const.h" #include "memory.h" #include "error.h" +#include "utils.h" +#include "fmt/format.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -62,27 +64,34 @@ FixQTB::FixQTB(LAMMPS *lmp, int narg, char **arg) : while (iarg < narg) { if (strcmp(arg[iarg],"temp") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qtb command"); - t_target = atof(arg[iarg+1]); if (t_target < 0.0) error->all(FLERR,"Fix qtb temp must be >= 0.0"); + t_target = force->numeric(FLERR,arg[iarg+1]); + if (t_target < 0.0) error->all(FLERR,"Fix qtb temp must be >= 0.0"); iarg += 2; } else if (strcmp(arg[iarg],"damp") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qtb command"); - t_period = atof(arg[iarg+1]); if (t_period <= 0.0) error->all(FLERR,"Fix qtb damp must be > 0.0"); fric_coef = 1/t_period; + t_period = force->numeric(FLERR,arg[iarg+1]); + if (t_period <= 0.0) error->all(FLERR,"Fix qtb damp must be > 0.0"); fric_coef = 1/t_period; iarg += 2; } else if (strcmp(arg[iarg],"seed") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qtb command"); - seed = atof(arg[iarg+1]); if (seed <= 0) error->all(FLERR,"Illegal fix qtb command"); + seed = force->inumeric(FLERR,arg[iarg+1]); + if (seed <= 0) error->all(FLERR,"Illegal fix qtb command"); iarg += 2; } else if (strcmp(arg[iarg],"f_max") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qtb command"); - f_max = atof(arg[iarg+1]); if (f_max <= 0) error->all(FLERR,"Illegal fix qtb command"); + f_max = force->numeric(FLERR,arg[iarg+1]); + if (f_max <= 0) error->all(FLERR,"Illegal fix qtb command"); iarg += 2; } else if (strcmp(arg[iarg],"N_f") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qtb command"); - N_f = atof(arg[iarg+1]); if (N_f <= 0) error->all(FLERR,"Illegal fix qtb command"); + N_f = force->inumeric(FLERR,arg[iarg+1]); + if (N_f <= 0) error->all(FLERR,"Illegal fix qtb command"); iarg += 2; } else error->all(FLERR,"Illegal fix qtb command"); } + maxexchange = 6*N_f+3; + // allocate qtb gfactor1 = NULL; gfactor3 = NULL; @@ -105,10 +114,6 @@ FixQTB::FixQTB(LAMMPS *lmp, int narg, char **arg) : // allocate random-arrays and fran grow_arrays(atom->nmax); atom->add_callback(0); - // memory->create(random_array_0,atom->nmax+300,2*N_f,"qtb:random_array_0"); - // memory->create(random_array_1,atom->nmax+300,2*N_f,"qtb:random_array_1"); - // memory->create(random_array_2,atom->nmax+300,2*N_f,"qtb:random_array_2"); - // memory->create(fran,atom->nmax+300,3,"qtb:fran"); // allocate omega_H and time_H memory->create(omega_H,2*N_f,"qtb:omega_H"); @@ -160,14 +165,17 @@ void FixQTB::init() //set up the h time step for updating the random force \delta{}h=\frac{\pi}{\Omega_{max}} if (int(1.0/(2*f_max*dtv)) == 0) { - if (comm->me == 0) printf ("Warning: Either f_max is too high or the time step is too big, setting f_max to be 1/timestep!\n"); + if (comm->me == 0) error->warning(FLERR,"Either f_max is too high or the time step " + "is too big, setting f_max to be 1/timestep!\n"); h_timestep=dtv; alpha=1; } else { alpha=int(1.0/(2*f_max*dtv)); h_timestep=alpha*dtv; } - if (comm->me == 0) printf ("The effective maximum frequncy is now %f inverse time unit with alpha value as %d!\n", 0.5/h_timestep, alpha); + if (comm->me == 0 && screen) + fmt::print(screen,"The effective maximum frequency is now {} inverse time unit " + "with alpha value as {}!\n", 0.5/h_timestep, alpha); // set force prefactors if (!atom->rmass) { @@ -404,11 +412,12 @@ void FixQTB::copy_arrays(int i, int j, int /*delflag*/) ------------------------------------------------------------------------- */ int FixQTB::pack_exchange(int i, double *buf) { - for (int m = 0; m < 2*N_f; m++) buf[m] = random_array_0[i][m]; - for (int m = 0; m < 2*N_f; m++) buf[m+2*N_f] = random_array_1[i][m]; - for (int m = 0; m < 2*N_f; m++) buf[m+4*N_f] = random_array_2[i][m]; - for (int m = 0; m < 3; m++) buf[m+6*N_f] = fran[i][m]; - return 6*N_f+3; + int n = 0; + for (int m = 0; m < 2*N_f; m++) buf[n++] = random_array_0[i][m]; + for (int m = 0; m < 2*N_f; m++) buf[n++] = random_array_1[i][m]; + for (int m = 0; m < 2*N_f; m++) buf[n++] = random_array_2[i][m]; + for (int m = 0; m < 3; m++) buf[n++] = fran[i][m]; + return n; } /* ---------------------------------------------------------------------- @@ -416,9 +425,10 @@ int FixQTB::pack_exchange(int i, double *buf) ------------------------------------------------------------------------- */ int FixQTB::unpack_exchange(int nlocal, double *buf) { - for (int m = 0; m < 2*N_f; m++) random_array_0[nlocal][m] = buf[m]; - for (int m = 0; m < 2*N_f; m++) random_array_1[nlocal][m] = buf[m+2*N_f]; - for (int m = 0; m < 2*N_f; m++) random_array_2[nlocal][m] = buf[m+4*N_f]; - for (int m = 0; m < 3; m++) fran[nlocal][m] = buf[m+6*N_f]; - return 6*N_f+3; + int n = 0; + for (int m = 0; m < 2*N_f; m++) random_array_0[nlocal][m] = buf[n++]; + for (int m = 0; m < 2*N_f; m++) random_array_1[nlocal][m] = buf[n++]; + for (int m = 0; m < 2*N_f; m++) random_array_2[nlocal][m] = buf[n++]; + for (int m = 0; m < 3; m++) fran[nlocal][m] = buf[n++]; + return n; } diff --git a/src/USER-QTB/fix_qtb.h b/src/USER-QTB/fix_qtb.h index 8fd14f5c5b..86a91bc9f6 100644 --- a/src/USER-QTB/fix_qtb.h +++ b/src/USER-QTB/fix_qtb.h @@ -47,22 +47,23 @@ class FixQTB : public Fix { private: // qtb parameters - int counter_mu; // counter l and mu - double t_period, fric_coef; // friction coefficient - int seed; // seed for the random number generator - double f_max; // frequency cutoff - int N_f; // number of frequency grid - double t_target; // target qtb temperature + int counter_mu; // counter l and mu + double t_period, fric_coef; // friction coefficient + int seed; // seed for the random number generator + double f_max; // frequency cutoff + int N_f; // number of frequency grid + double t_target; // target qtb temperature char *id_temp; class Compute *temperature; - double h_timestep; // time step to update the random forces - int alpha; // number of time steps to update the random forces - class RanMars *random; // random number generator - double *gfactor1,*gfactor3; // factors of frictions and random forces - double *omega_H,*time_H; // H gives the desired power spectrum - double **random_array_0, **random_array_1, **random_array_2; // random number arrays give independence between atoms and directions + double h_timestep; // time step to update the random forces + int alpha; // number of time steps to update the random forces + class RanMars *random; // random number generator + double *gfactor1,*gfactor3; // factors of frictions and random forces + double *omega_H,*time_H; // H gives the desired power spectrum + // random number arrays give independence between atoms and directions + double **random_array_0, **random_array_1, **random_array_2; int nlevels_respa; - double **fran, fsum[3], fsumall[3]; // random forces and their sums + double **fran, fsum[3], fsumall[3]; // random forces and their sums }; } diff --git a/src/USER-REAXC/reaxc_ffield.cpp b/src/USER-REAXC/reaxc_ffield.cpp index c8e097eb1c..b865b48ad3 100644 --- a/src/USER-REAXC/reaxc_ffield.cpp +++ b/src/USER-REAXC/reaxc_ffield.cpp @@ -337,6 +337,8 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax, j = atoi(tmp[0]) - 1; k = atoi(tmp[1]) - 1; + if ((c < 10) || (j < 0) || (k < 0)) + control->error_ptr->all(FLERR, "Inconsistent force field file"); if (j < reax->num_atom_types && k < reax->num_atom_types) { @@ -361,6 +363,8 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax, /* line 2 */ fgets(s,MAX_LINE,fp); c=Tokenize(s,&tmp); + if (c < 8) + control->error_ptr->all(FLERR, "Inconsistent force field file"); val = atof(tmp[0]); reax->tbp[j][k].p_be2 = val; reax->tbp[k][j].p_be2 = val; @@ -491,6 +495,9 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax, j = atoi(tmp[0]) - 1; k = atoi(tmp[1]) - 1; + if ((c < (lgflag ? 9 : 8)) || (j < 0) || (k < 0)) + control->error_ptr->all(FLERR, "Inconsistent force field file"); + if (j < reax->num_atom_types && k < reax->num_atom_types) { val = atof(tmp[2]); if (val > 0.0) { @@ -528,10 +535,12 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax, reax->tbp[k][j].r_pp = val; } - val = atof(tmp[8]); - if (val >= 0.0) { - reax->tbp[j][k].lgcij = val; - reax->tbp[k][j].lgcij = val; + if (lgflag) { + val = atof(tmp[8]); + if (val >= 0.0) { + reax->tbp[j][k].lgcij = val; + reax->tbp[k][j].lgcij = val; + } } } } @@ -552,6 +561,8 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax, j = atoi(tmp[0]) - 1; k = atoi(tmp[1]) - 1; m = atoi(tmp[2]) - 1; + if ((c < 10) || (j < 0) || (k < 0) || (m < 0)) + control->error_ptr->all(FLERR, "Inconsistent force field file"); if (j < reax->num_atom_types && k < reax->num_atom_types && m < reax->num_atom_types) { @@ -611,6 +622,8 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax, k = atoi(tmp[1]) - 1; m = atoi(tmp[2]) - 1; n = atoi(tmp[3]) - 1; + if ((c < 9) || (k < 0) || (m < 0)) + control->error_ptr->all(FLERR, "Inconsistent force field file"); if (j >= 0 && n >= 0) { // this means the entry is not in compact form if (j < reax->num_atom_types && k < reax->num_atom_types && @@ -667,8 +680,6 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax, } } - - /* next line is number of hydrogen bond params and some comments */ fgets( s, MAX_LINE, fp ); c = Tokenize( s, &tmp ); @@ -686,7 +697,8 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax, j = atoi(tmp[0]) - 1; k = atoi(tmp[1]) - 1; m = atoi(tmp[2]) - 1; - + if ((c < 7) || (j < 0) || (k < 0) || (m < 0)) + control->error_ptr->all(FLERR, "Inconsistent force field file"); if (j < reax->num_atom_types && m < reax->num_atom_types) { val = atof(tmp[3]); diff --git a/src/atom.cpp b/src/atom.cpp index f47906f203..f6379213ea 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -48,8 +48,6 @@ using namespace MathConst; #define DELTA_PERATOM 64 #define EPSILON 1.0e-6 -enum{DOUBLE,INT,BIGINT}; - /* ---------------------------------------------------------------------- */ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) diff --git a/src/atom.h b/src/atom.h index fca682aff6..e20ee0596d 100644 --- a/src/atom.h +++ b/src/atom.h @@ -24,6 +24,7 @@ class Atom : protected Pointers { public: char *atom_style; class AtomVec *avec; + enum{DOUBLE,INT,BIGINT}; // atom counts diff --git a/src/atom_vec.cpp b/src/atom_vec.cpp index b9a8d6b09d..8cd7db4fd4 100644 --- a/src/atom_vec.cpp +++ b/src/atom_vec.cpp @@ -33,8 +33,6 @@ using namespace MathConst; #define DELTA 16384 #define DELTA_BONUS 8192 -enum{DOUBLE,INT,BIGINT}; - /* ---------------------------------------------------------------------- */ AtomVec::AtomVec(LAMMPS *lmp) : Pointers(lmp) @@ -214,7 +212,7 @@ void AtomVec::grow(int n) datatype = mgrow.datatype[i]; cols = mgrow.cols[i]; const int nthreads = threads[i] ? comm->nthreads : 1; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) memory->grow(*((double **) pdata),nmax*nthreads,"atom:dvec"); else if (cols > 0) @@ -223,7 +221,7 @@ void AtomVec::grow(int n) maxcols = *(mgrow.maxcols[i]); memory->grow(*((double ***) pdata),nmax*nthreads,maxcols,"atom:darray"); } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) memory->grow(*((int **) pdata),nmax*nthreads,"atom:ivec"); else if (cols > 0) @@ -232,14 +230,14 @@ void AtomVec::grow(int n) maxcols = *(mgrow.maxcols[i]); memory->grow(*((int ***) pdata),nmax*nthreads,maxcols,"atom:iarray"); } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) memory->grow(*((bigint **) pdata),nmax*nthreads,"atom:bvec"); else if (cols > 0) memory->grow(*((bigint ***) pdata),nmax*nthreads,cols,"atom:barray"); else { maxcols = *(mgrow.maxcols[i]); - memory->grow(*((int ***) pdata),nmax*nthreads,maxcols,"atom:barray"); + memory->grow(*((bigint ***) pdata),nmax*nthreads,maxcols,"atom:barray"); } } } @@ -275,7 +273,7 @@ void AtomVec::copy(int i, int j, int delflag) pdata = mcopy.pdata[n]; datatype = mcopy.datatype[n]; cols = mcopy.cols[n]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); vec[j] = vec[i]; @@ -292,7 +290,7 @@ void AtomVec::copy(int i, int j, int delflag) for (m = 0; m < ncols; m++) array[j][m] = array[i][m]; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); vec[j] = vec[i]; @@ -309,7 +307,7 @@ void AtomVec::copy(int i, int j, int delflag) for (m = 0; m < ncols; m++) array[j][m] = array[i][m]; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); vec[j] = vec[i]; @@ -377,7 +375,7 @@ int AtomVec::pack_comm(int n, int *list, double *buf, pdata = mcomm.pdata[nn]; datatype = mcomm.datatype[nn]; cols = mcomm.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); for (i = 0; i < n; i++) { @@ -392,7 +390,7 @@ int AtomVec::pack_comm(int n, int *list, double *buf, buf[m++] = array[j][mm]; } } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); for (i = 0; i < n; i++) { @@ -407,7 +405,7 @@ int AtomVec::pack_comm(int n, int *list, double *buf, buf[m++] = ubuf(array[j][mm]).d; } } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); for (i = 0; i < n; i++) { @@ -498,7 +496,7 @@ int AtomVec::pack_comm_vel(int n, int *list, double *buf, pdata = mcomm_vel.pdata[nn]; datatype = mcomm_vel.datatype[nn]; cols = mcomm_vel.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); for (i = 0; i < n; i++) { @@ -513,7 +511,7 @@ int AtomVec::pack_comm_vel(int n, int *list, double *buf, buf[m++] = array[j][mm]; } } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); for (i = 0; i < n; i++) { @@ -528,7 +526,7 @@ int AtomVec::pack_comm_vel(int n, int *list, double *buf, buf[m++] = ubuf(array[j][mm]).d; } } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); for (i = 0; i < n; i++) { @@ -572,7 +570,7 @@ void AtomVec::unpack_comm(int n, int first, double *buf) pdata = mcomm.pdata[nn]; datatype = mcomm.datatype[nn]; cols = mcomm.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); for (i = first; i < last; i++) @@ -583,7 +581,7 @@ void AtomVec::unpack_comm(int n, int first, double *buf) for (mm = 0; mm < cols; mm++) array[i][mm] = buf[m++]; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); for (i = first; i < last; i++) @@ -594,7 +592,7 @@ void AtomVec::unpack_comm(int n, int first, double *buf) for (mm = 0; mm < cols; mm++) array[i][mm] = (int) ubuf(buf[m++]).i; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); for (i = first; i < last; i++) @@ -635,7 +633,7 @@ void AtomVec::unpack_comm_vel(int n, int first, double *buf) pdata = mcomm_vel.pdata[nn]; datatype = mcomm_vel.datatype[nn]; cols = mcomm_vel.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); for (i = first; i < last; i++) @@ -646,18 +644,18 @@ void AtomVec::unpack_comm_vel(int n, int first, double *buf) for (mm = 0; mm < cols; mm++) array[i][mm] = buf[m++]; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); for (i = first; i < last; i++) - vec[i] = ubuf(buf[m++]).i; + vec[i] = (int) ubuf(buf[m++]).i; } else { int **array = *((int ***) pdata); for (i = first; i < last; i++) for (mm = 0; mm < cols; mm++) array[i][mm] = (int) ubuf(buf[m++]).i; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); for (i = first; i < last; i++) @@ -695,7 +693,7 @@ int AtomVec::pack_reverse(int n, int first, double *buf) pdata = mreverse.pdata[nn]; datatype = mreverse.datatype[nn]; cols = mreverse.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); for (i = first; i < last; i++) { @@ -708,7 +706,7 @@ int AtomVec::pack_reverse(int n, int first, double *buf) buf[m++] = array[i][mm]; } } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); for (i = first; i < last; i++) { @@ -721,7 +719,7 @@ int AtomVec::pack_reverse(int n, int first, double *buf) buf[m++] = ubuf(array[i][mm]).d; } } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); for (i = first; i < last; i++) { @@ -761,7 +759,7 @@ void AtomVec::unpack_reverse(int n, int *list, double *buf) pdata = mreverse.pdata[nn]; datatype = mreverse.datatype[nn]; cols = mreverse.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); for (i = 0; i < n; i++) { @@ -776,34 +774,34 @@ void AtomVec::unpack_reverse(int n, int *list, double *buf) array[j][mm] += buf[m++]; } } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); for (i = 0; i < n; i++) { j = list[i]; - vec[j] += buf[m++]; + vec[j] += (int) ubuf(buf[m++]).i; } } else { int **array = *((int ***) pdata); for (i = 0; i < n; i++) { j = list[i]; for (mm = 0; mm < cols; mm++) - array[j][mm] += buf[m++]; + array[j][mm] += (int) ubuf(buf[m++]).i; } } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); for (i = 0; i < n; i++) { j = list[i]; - vec[j] += buf[m++]; + vec[j] += (bigint) ubuf(buf[m++]).i; } } else { bigint **array = *((bigint ***) pdata); for (i = 0; i < n; i++) { j = list[i]; for (mm = 0; mm < cols; mm++) - array[j][mm] += buf[m++]; + array[j][mm] += (bigint) ubuf(buf[m++]).i; } } } @@ -856,7 +854,7 @@ int AtomVec::pack_border(int n, int *list, double *buf, int pbc_flag, int *pbc) pdata = mborder.pdata[nn]; datatype = mborder.datatype[nn]; cols = mborder.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); for (i = 0; i < n; i++) { @@ -871,7 +869,7 @@ int AtomVec::pack_border(int n, int *list, double *buf, int pbc_flag, int *pbc) buf[m++] = array[j][mm]; } } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); for (i = 0; i < n; i++) { @@ -886,7 +884,7 @@ int AtomVec::pack_border(int n, int *list, double *buf, int pbc_flag, int *pbc) buf[m++] = ubuf(array[j][mm]).d; } } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); for (i = 0; i < n; i++) { @@ -990,7 +988,7 @@ int AtomVec::pack_border_vel(int n, int *list, double *buf, pdata = mborder_vel.pdata[nn]; datatype = mborder_vel.datatype[nn]; cols = mborder_vel.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); for (i = 0; i < n; i++) { @@ -1005,7 +1003,7 @@ int AtomVec::pack_border_vel(int n, int *list, double *buf, buf[m++] = array[j][mm]; } } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); for (i = 0; i < n; i++) { @@ -1020,7 +1018,7 @@ int AtomVec::pack_border_vel(int n, int *list, double *buf, buf[m++] = ubuf(array[j][mm]).d; } } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); for (i = 0; i < n; i++) { @@ -1072,7 +1070,7 @@ void AtomVec::unpack_border(int n, int first, double *buf) pdata = mborder.pdata[nn]; datatype = mborder.datatype[nn]; cols = mborder.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); for (i = first; i < last; i++) @@ -1083,7 +1081,7 @@ void AtomVec::unpack_border(int n, int first, double *buf) for (mm = 0; mm < cols; mm++) array[i][mm] = buf[m++]; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); for (i = first; i < last; i++) @@ -1094,7 +1092,7 @@ void AtomVec::unpack_border(int n, int first, double *buf) for (mm = 0; mm < cols; mm++) array[i][mm] = (int) ubuf(buf[m++]).i; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); for (i = first; i < last; i++) @@ -1144,7 +1142,7 @@ void AtomVec::unpack_border_vel(int n, int first, double *buf) pdata = mborder_vel.pdata[nn]; datatype = mborder_vel.datatype[nn]; cols = mborder_vel.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); for (i = first; i < last; i++) @@ -1155,18 +1153,18 @@ void AtomVec::unpack_border_vel(int n, int first, double *buf) for (mm = 0; mm < cols; mm++) array[i][mm] = buf[m++]; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); for (i = first; i < last; i++) - vec[i] = ubuf(buf[m++]).i; + vec[i] = (int) ubuf(buf[m++]).i; } else { int **array = *((int ***) pdata); for (i = first; i < last; i++) for (mm = 0; mm < cols; mm++) array[i][mm] = (int) ubuf(buf[m++]).i; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); for (i = first; i < last; i++) @@ -1216,7 +1214,7 @@ int AtomVec::pack_exchange(int i, double *buf) pdata = mexchange.pdata[nn]; datatype = mexchange.datatype[nn]; cols = mexchange.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); buf[m++] = vec[i]; @@ -1233,7 +1231,7 @@ int AtomVec::pack_exchange(int i, double *buf) for (mm = 0; mm < ncols; mm++) buf[m++] = array[i][mm]; } - } if (datatype == INT) { + } if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); buf[m++] = ubuf(vec[i]).d; @@ -1250,7 +1248,7 @@ int AtomVec::pack_exchange(int i, double *buf) for (mm = 0; mm < ncols; mm++) buf[m++] = ubuf(array[i][mm]).d; } - } if (datatype == BIGINT) { + } if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); buf[m++] = ubuf(vec[i]).d; @@ -1308,7 +1306,7 @@ int AtomVec::unpack_exchange(double *buf) pdata = mexchange.pdata[nn]; datatype = mexchange.datatype[nn]; cols = mexchange.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); vec[nlocal] = buf[m++]; @@ -1325,10 +1323,10 @@ int AtomVec::unpack_exchange(double *buf) for (mm = 0; mm < ncols; mm++) array[nlocal][mm] = buf[m++]; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); - vec[nlocal] = ubuf(buf[m++]).i; + vec[nlocal] = (int) ubuf(buf[m++]).i; } else if (cols > 0) { int **array = *((int ***) pdata); for (mm = 0; mm < cols; mm++) @@ -1342,7 +1340,7 @@ int AtomVec::unpack_exchange(double *buf) for (mm = 0; mm < ncols; mm++) array[nlocal][mm] = (int) ubuf(buf[m++]).i; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); vec[nlocal] = (bigint) ubuf(buf[m++]).i; @@ -1450,7 +1448,7 @@ int AtomVec::pack_restart(int i, double *buf) pdata = mrestart.pdata[nn]; datatype = mrestart.datatype[nn]; cols = mrestart.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); buf[m++] = vec[i]; @@ -1467,7 +1465,7 @@ int AtomVec::pack_restart(int i, double *buf) for (mm = 0; mm < ncols; mm++) buf[m++] = array[i][mm]; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); buf[m++] = ubuf(vec[i]).d; @@ -1484,7 +1482,7 @@ int AtomVec::pack_restart(int i, double *buf) for (mm = 0; mm < ncols; mm++) buf[m++] = ubuf(array[i][mm]).d; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); buf[m++] = ubuf(vec[i]).d; @@ -1551,7 +1549,7 @@ int AtomVec::unpack_restart(double *buf) pdata = mrestart.pdata[nn]; datatype = mrestart.datatype[nn]; cols = mrestart.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); vec[nlocal] = buf[m++]; @@ -1568,10 +1566,10 @@ int AtomVec::unpack_restart(double *buf) for (mm = 0; mm < ncols; mm++) array[nlocal][mm] = buf[m++]; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); - vec[nlocal] = ubuf(buf[m++]).i; + vec[nlocal] = (int) ubuf(buf[m++]).i; } else if (cols > 0) { int **array = *((int ***) pdata); for (mm = 0; mm < cols; mm++) @@ -1585,7 +1583,7 @@ int AtomVec::unpack_restart(double *buf) for (mm = 0; mm < ncols; mm++) array[nlocal][mm] = (int) ubuf(buf[m++]).i; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); vec[nlocal] = (bigint) ubuf(buf[m++]).i; @@ -1594,7 +1592,7 @@ int AtomVec::unpack_restart(double *buf) for (mm = 0; mm < cols; mm++) array[nlocal][mm] = (bigint) ubuf(buf[m++]).i; } else { - int **array = *((int ***) pdata); + bigint **array = *((bigint ***) pdata); collength = mexchange.collength[nn]; plength = mexchange.plength[nn]; if (collength) ncols = (*((int ***) plength))[nlocal][collength-1]; @@ -1654,7 +1652,7 @@ void AtomVec::create_atom(int itype, double *coord) pdata = mcreate.pdata[n]; datatype = mcreate.datatype[n]; cols = mcreate.cols[n]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); vec[nlocal] = 0.0; @@ -1663,7 +1661,7 @@ void AtomVec::create_atom(int itype, double *coord) for (m = 0; m < cols; m++) array[nlocal][m] = 0.0; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); vec[nlocal] = 0; @@ -1672,7 +1670,7 @@ void AtomVec::create_atom(int itype, double *coord) for (m = 0; m < cols; m++) array[nlocal][m] = 0; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); vec[nlocal] = 0; @@ -1718,7 +1716,7 @@ void AtomVec::data_atom(double *coord, imageint imagetmp, char **values) pdata = mdata_atom.pdata[n]; datatype = mdata_atom.datatype[n]; cols = mdata_atom.cols[n]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); vec[nlocal] = utils::numeric(FLERR,values[ivalue++],true,lmp); @@ -1731,7 +1729,7 @@ void AtomVec::data_atom(double *coord, imageint imagetmp, char **values) for (m = 0; m < cols; m++) array[nlocal][m] = utils::numeric(FLERR,values[ivalue++],true,lmp); } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); vec[nlocal] = utils::inumeric(FLERR,values[ivalue++],true,lmp); @@ -1740,7 +1738,7 @@ void AtomVec::data_atom(double *coord, imageint imagetmp, char **values) for (m = 0; m < cols; m++) array[nlocal][m] = utils::inumeric(FLERR,values[ivalue++],true,lmp); } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); vec[nlocal] = utils::bnumeric(FLERR,values[ivalue++],true,lmp); @@ -1788,7 +1786,7 @@ void AtomVec::pack_data(double **buf) pdata = mdata_atom.pdata[n]; datatype = mdata_atom.datatype[n]; cols = mdata_atom.cols[n]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); buf[i][j++] = vec[i]; @@ -1797,7 +1795,7 @@ void AtomVec::pack_data(double **buf) for (m = 0; m < cols; m++) buf[i][j++] = array[i][m]; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); buf[i][j++] = ubuf(vec[i]).d; @@ -1806,7 +1804,7 @@ void AtomVec::pack_data(double **buf) for (m = 0; m < cols; m++) buf[i][j++] = ubuf(array[i][m]).d; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); buf[i][j++] = ubuf(vec[i]).d; @@ -1844,21 +1842,21 @@ void AtomVec::write_data(FILE *fp, int n, double **buf) for (nn = 1; nn < ndata_atom; nn++) { datatype = mdata_atom.datatype[nn]; cols = mdata_atom.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { fprintf(fp," %-1.16e",buf[i][j++]); } else { for (m = 0; m < cols; m++) fprintf(fp," %-1.16e",buf[i][j++]); } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { fprintf(fp," %d",(int) ubuf(buf[i][j++]).i); } else { for (m = 0; m < cols; m++) fprintf(fp," %d",(int) ubuf(buf[i][j++]).i); } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { fmt::print(fp," {}",(bigint) ubuf(buf[i][j++]).i); } else { @@ -1895,7 +1893,7 @@ void AtomVec::data_vel(int ilocal, char **values) pdata = mdata_vel.pdata[n]; datatype = mdata_vel.datatype[n]; cols = mdata_vel.cols[n]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); vec[ilocal] = utils::numeric(FLERR,values[ivalue++],true,lmp); @@ -1904,7 +1902,7 @@ void AtomVec::data_vel(int ilocal, char **values) for (m = 0; m < cols; m++) array[ilocal][m] = utils::numeric(FLERR,values[ivalue++],true,lmp); } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); vec[ilocal] = utils::inumeric(FLERR,values[ivalue++],true,lmp); @@ -1913,7 +1911,7 @@ void AtomVec::data_vel(int ilocal, char **values) for (m = 0; m < cols; m++) array[ilocal][m] = utils::inumeric(FLERR,values[ivalue++],true,lmp); } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); vec[ilocal] = utils::bnumeric(FLERR,values[ivalue++],true,lmp); @@ -1944,7 +1942,7 @@ void AtomVec::pack_vel(double **buf) pdata = mdata_vel.pdata[n]; datatype = mdata_vel.datatype[n]; cols = mdata_vel.cols[n]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { double *vec = *((double **) pdata); buf[i][j++] = vec[i]; @@ -1953,7 +1951,7 @@ void AtomVec::pack_vel(double **buf) for (m = 0; m < cols; m++) buf[i][j++] = array[i][m]; } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { int *vec = *((int **) pdata); buf[i][j++] = ubuf(vec[i]).d; @@ -1962,7 +1960,7 @@ void AtomVec::pack_vel(double **buf) for (m = 0; m < cols; m++) buf[i][j++] = ubuf(array[i][m]).d; } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bigint *vec = *((bigint **) pdata); buf[i][j++] = ubuf(vec[i]).d; @@ -1992,21 +1990,21 @@ void AtomVec::write_vel(FILE *fp, int n, double **buf) for (nn = 1; nn < ndata_vel; nn++) { datatype = mdata_vel.datatype[nn]; cols = mdata_vel.cols[nn]; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { fprintf(fp," %-1.16e",buf[i][j++]); } else { for (m = 0; m < cols; m++) fprintf(fp," %-1.16e",buf[i][j++]); } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { fprintf(fp," %d",(int) ubuf(buf[i][j++]).i); } else { for (m = 0; m < cols; m++) fprintf(fp," %d",(int) ubuf(buf[i][j++]).i); } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { fmt::print(fp," {}",(bigint) ubuf(buf[i][j++]).i); } else { @@ -2290,7 +2288,7 @@ bigint AtomVec::memory_usage() datatype = mgrow.datatype[i]; cols = mgrow.cols[i]; const int nthreads = threads[i] ? comm->nthreads : 1; - if (datatype == DOUBLE) { + if (datatype == Atom::DOUBLE) { if (cols == 0) { bytes += memory->usage(*((double **) pdata),nmax*nthreads); } else if (cols > 0) { @@ -2299,7 +2297,7 @@ bigint AtomVec::memory_usage() maxcols = *(mgrow.maxcols[i]); bytes += memory->usage(*((double ***) pdata),nmax*nthreads,maxcols); } - } else if (datatype == INT) { + } else if (datatype == Atom::INT) { if (cols == 0) { bytes += memory->usage(*((int **) pdata),nmax*nthreads); } else if (cols > 0) { @@ -2308,7 +2306,7 @@ bigint AtomVec::memory_usage() maxcols = *(mgrow.maxcols[i]); bytes += memory->usage(*((int ***) pdata),nmax*nthreads,maxcols); } - } else if (datatype == BIGINT) { + } else if (datatype == Atom::BIGINT) { if (cols == 0) { bytes += memory->usage(*((bigint **) pdata),nmax*nthreads); } else if (cols > 0) { diff --git a/src/comm.cpp b/src/comm.cpp index 3df3cc4d44..9102e2991d 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -258,10 +258,8 @@ void Comm::init_exchange() int onefix; maxexchange_fix = 0; - for (int i = 0; i < nfix; i++) { - onefix = fix[i]->maxexchange; - maxexchange_fix = MAX(maxexchange_fix,onefix); - } + for (int i = 0; i < nfix; i++) + maxexchange_fix += fix[i]->maxexchange; maxexchange = maxexchange_atom + maxexchange_fix; bufextra = maxexchange + BUFEXTRA; diff --git a/src/force.cpp b/src/force.cpp index 2827c43268..840d7e31ca 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -1019,7 +1019,8 @@ FILE *Force::open_potential(const char *name) date = utils::get_potential_date(filepath, "potential"); if(!date.empty()) { - utils::logmesg(lmp, fmt::format("Reading potential file {} with DATE: {}", name, date)); + utils::logmesg(lmp, fmt::format("Reading potential file {} " + "with DATE: {}\n", name, date)); } return fopen(filepath.c_str(), "r"); } diff --git a/tools/singularity/centos7.def b/tools/singularity/centos7.def index 7803f3f386..5369c6977b 100644 --- a/tools/singularity/centos7.def +++ b/tools/singularity/centos7.def @@ -11,9 +11,43 @@ From: centos:7 hdf5-devel python36-virtualenv python36-pip python-pip \ netcdf-devel netcdf-cxx-devel netcdf-mpich-devel netcdf-openmpi-devel \ python-virtualenv fftw-devel voro++-devel eigen3-devel gsl-devel openblas-devel enchant \ - blas-devel lapack-devel libyaml-devel + blas-devel lapack-devel libyaml-devel openkim-models kim-api-devel yum clean all + # we need to reset any module variables + # inherited from the host. + unset __LMOD_REF_COUNT__LMFILES_ + unset __LMOD_REF_COUNT_PATH + unset __LMOD_REF_COUNT_LD_LIBRARY_PATH + unset __LMOD_REF_COUNT_MANPATH + unset __LMOD_REF_COUNT_MODULEPATH + unset __LMOD_REF_COUNT_LOADEDMODULES + unset _LMFILES_ + unset MODULEPATH + unset MODULESHOME + unset MODULEPATH_ROOT + unset LOADEDMODULES + unset LMOD_SYSTEM_DEFAULT_MODULES + + # load MPI by default + . /etc/profile + module load mpi + + # manually install Plumed + mkdir plumed + cd plumed + version=2.6.0 + curl -L -o plumed.tar.gz https://github.com/plumed/plumed2/releases/download/v${version}/plumed-src-${version}.tgz + tar -xzf plumed.tar.gz + cd plumed-${version} + ./configure --disable-doc --prefix=/usr + make + make install + # fix up installation for CentOS and Fedora + mv -v /usr/lib/pkgconfig/plumed* /usr/share/pkgconfig/ + cd ../../ + rm -rvf plumed + # set custom prompt indicating the container name CUSTOM_PROMPT_ENV=/.singularity.d/env/99-zz_custom_prompt.sh cat >$CUSTOM_PROMPT_ENV <$CUSTOM_PROMPT_ENV <$CUSTOM_PROMPT_ENV <$CUSTOM_PROMPT_ENV <$CUSTOM_PROMPT_ENV <