Merge pull request #2132 from akohlmey/collected-bugfixes-and-updates

Collected bugfixes and updates
This commit is contained in:
Axel Kohlmeyer
2020-06-08 19:55:34 -04:00
committed by GitHub
23 changed files with 394 additions and 253 deletions

View File

@ -94,6 +94,7 @@ An alphabetic list of all general LAMMPS commands.
* :doc:`package <package>`
* :doc:`pair_coeff <pair_coeff>`
* :doc:`pair_modify <pair_modify>`
* :doc:`pair_style <pair_style>`
* :doc:`pair_write <pair_write>`
* :doc:`partition <partition>`
* :doc:`prd <prd>`

View File

@ -57,6 +57,8 @@
------------------------------------------------------------------------- */
#include "kim_init.h"
#include "fix_store_kim.h"
#include "kim_units.h"
#include <cstring>
#include <string>
#include <sstream>
@ -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"

View File

@ -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);

View File

@ -21,6 +21,7 @@
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <string>
#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;
}

View File

@ -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;
}

View File

@ -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
};
}

View File

@ -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]);

View File

@ -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)

View File

@ -24,6 +24,7 @@ class Atom : protected Pointers {
public:
char *atom_style;
class AtomVec *avec;
enum{DOUBLE,INT,BIGINT};
// atom counts

View File

@ -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) {

View File

@ -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;

View File

@ -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");
}

View File

@ -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 <<EOF

View File

@ -16,9 +16,43 @@ From: centos:8
texlive-latex-bin texlive-lualatex-math texlive-fncychap texlive-tabulary \
texlive-framed texlive-wrapfig texlive-upquote texlive-capt-of \
texlive-needspace texlive-titlesec texlive-anysize texlive-dvipng \
blas-devel lapack-devel libyaml-devel
blas-devel lapack-devel libyaml-devel openkim-models kim-api-devel
dnf 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/lib64/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 <<EOF

View File

@ -35,9 +35,29 @@ From: fedora:32
texlive-latex-bin texlive-lualatex-math texlive-fncychap texlive-tabulary \
texlive-framed texlive-wrapfig texlive-upquote texlive-capt-of \
texlive-needspace texlive-titlesec texlive-anysize texlive-dvipng \
blas-devel lapack-devel libyaml-devel
blas-devel lapack-devel libyaml-devel openkim-models kim-api-devel
dnf clean all
# enable Lmod and load MPI
source /usr/share/lmod/lmod/init/profile
module purge
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/lib64/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 <<EOF

View File

@ -4,6 +4,9 @@ From: ubuntu:18.04
%post
export DEBIAN_FRONTEND=noninteractive
apt-get update
apt-get install --no-install-recommends -y software-properties-common
add-apt-repository ppa:openkim/latest
apt-get update
apt-get upgrade --no-install-recommends -y
apt-get install --no-install-recommends -y \
bc \
@ -66,11 +69,28 @@ From: ubuntu:18.04
wget \
xxd \
valgrind \
gdb
gdb \
libkim-api-dev \
openkim-models
# clean cache
rm -rf /var/lib/apt/lists/*
# 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 <<EOF

View File

@ -4,6 +4,9 @@ From: ubuntu:20.04
%post
export DEBIAN_FRONTEND=noninteractive
apt-get update
apt-get install --no-install-recommends -y software-properties-common
add-apt-repository ppa:openkim/latest
apt-get update
apt-get upgrade --no-install-recommends -y
apt-get install --no-install-recommends -y \
bc \
@ -62,11 +65,28 @@ From: ubuntu:20.04
wget \
xxd \
valgrind \
gdb
gdb \
libkim-api-dev \
openkim-models
# clean cache
rm -rf /var/lib/apt/lists/*
# 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 <<EOF
@ -75,10 +95,10 @@ PS1="[ubuntu20.04:\u@\h] \W> "
EOF
chmod 755 $CUSTOM_PROMPT_ENV
%environment
LC_ALL=C
export LC_ALL
export PATH=/usr/lib/ccache:$PATH
# restrict OpenMPI to shared memory comm by default
OMPI_MCA_btl="tcp,self"
# do not warn about unused components as this messes up testing

View File

@ -1,7 +1,7 @@
---
lammps_version: 5 May 2020
date_generated: Sat May 30 17:49:16 202
epsilon: 2.5e-13
epsilon: 5e-13
prerequisites: ! |
atom full
bond gromos

View File

@ -1,7 +1,7 @@
---
lammps_version: 5 May 2020
date_generated: Sat May 30 17:49:16 202
epsilon: 2.5e-13
epsilon: 5e-13
prerequisites: ! |
atom full
bond nonlinear

View File

@ -1,7 +1,7 @@
---
lammps_version: 5 May 2020
date_generated: Sun May 31 08:39:49 202
epsilon: 5e-14
epsilon: 1e-13
prerequisites: ! |
atom full
pair coul/cut

View File

@ -1,7 +1,7 @@
---
lammps_version: 5 May 2020
date_generated: Sun May 31 08:50:20 202
epsilon: 1e-13
epsilon: 2e-13
prerequisites: ! |
atom full
pair coul/dsf

View File

@ -1,7 +1,7 @@
---
lammps_version: 5 May 2020
date_generated: Sun May 31 09:23:33 202
epsilon: 1e-13
epsilon: 5e-13
prerequisites: ! |
atom full
pair coul/wolf

View File

@ -1,7 +1,7 @@
---
lammps_version: 5 May 2020
date_generated: Sun May 31 07:05:48 202
epsilon: 5e-14
epsilon: 2e-13
prerequisites: ! |
atom full
pair lj/charmm/coul/charmm/implicit