ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'develop' of github.com:lammps/lammps into skip_trim

Browse Source
This commit is contained in:
Stan Gerald Moore
2023-11-17 17:49:27 -07:00
parent be19b5c210 9a432341f9
commit d1f3411db9
5 changed files with 277 additions and 93 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
doc/src/Howto_body.rst
View File

@ -170,9 +170,9 @@ with this body style to compute body/body and body/non-body interactions.
The *rounded/polygon* body style represents body particles as a 2d
polygon with a variable number of N vertices. This style can only be
used for 2d models; see the :doc:`boundary <boundary>` command. See the
"pair_style body/rounded/polygon" page for a diagram of two
squares with rounded circles at the vertices. Special cases for N = 1
(circle) and N = 2 (rod with rounded ends) can also be specified.
:doc:`pair_style body/rounded/polygon <pair_body_rounded_polygon>` page for
a diagram of two squares with rounded circles at the vertices. Special cases
for N = 1 (circle) and N = 2 (rod with rounded ends) can also be specified.
One use of this body style is for 2d discrete element models, as
described in :ref:`Fraige <body-Fraige>`.

16
doc/src/fix_pimd.rst
View File

@ -31,7 +31,7 @@ Syntax
.. parsed-literal::
*keywords* = *method* or *integrator* or *ensemble* or *fmmode* or *fmass* or *scale* or *temp* or *thermostat* or *tau* or *iso* or *aniso* or *barostat* or *taup* or *fixcom* or *lj*
*method* value = *nmpimd*
*method* value = *nmpimd* (default) or *pimd*
*integrator* value = *obabo* or *baoab*
*fmmode* value = *physical* or *normal*
*fmass* value = scaling factor on mass
@ -137,8 +137,6 @@ normal-mode PIMD. A value of *cmd* is for centroid molecular dynamics
the real particle.
.. note::
Fix pimd/langevin only supports *method* value *nmpimd*. This should be enough
for most PIMD applications for quantum thermodynamics purpose.
Motion of the centroid can be effectively uncoupled from the other
normal modes by scaling the fictitious masses to achieve a partial
@ -151,6 +149,16 @@ normal-mode PIMD. A value of *cmd* is for centroid molecular dynamics
only the k > 0 modes are thermostatted, not the centroid degrees of
freedom.
.. versionadded:: TBD
Mode *pimd* added to fix pimd/langevin.
Fix pimd/langevin supports the *method* values *nmpimd* and *pimd*. The default value is *nmpimd*.
If *method* is *nmpimd*, the normal mode representation is used to integrate the equations of motion.
The exact solution of harmonic oscillator is used to propagate the free ring polymer part of the Hamiltonian.
If *method* is *pimd*, the Cartesian representation is used to integrate the equations of motion.
The harmonic force is added to the total force of the system, and the numerical integrator is used to propagate the Hamiltonian.
The keyword *integrator* specifies the Trotter splitting method used by *fix pimd/langevin*.
See :ref:`(Liu) <Liu>` for a discussion on the OBABO and BAOAB splitting schemes. Typically
either of the two should work fine.
@ -207,6 +215,7 @@ The keyword *thermostat* reads *style* and *seed* of thermostat for fix style *p
be *PILE_L* (path integral Langevin equation local thermostat, as described in :ref:`Ceriotti <Ceriotti2>`), and *seed* should a positive integer number, which serves as the seed of the pseudo random number generator.
.. note::
The fix style *pimd/langevin* uses the stochastic PILE_L thermostat to control temperature. This thermostat works on the normal modes
of the ring polymer. The *tau* parameter controls the centroid mode, and the *scale* parameter controls the non-centroid modes.
@ -269,6 +278,7 @@ related tasks for each of the partitions, e.g.
read_restart system_${ibead}.restart2
.. note::
Fix *pimd/langevin* dumps the Cartesian coordinates, but dumps the velocities and
forces in the normal mode representation. If the Cartesian velocities and forces are
needed, it is easy to perform the transformation when doing post-processing.

2
doc/src/pair_body_rounded_polyhedron.rst
View File

@ -40,7 +40,7 @@ rounded/polyhedron particles.
This pairwise interaction between the rounded polyhedra is described
in :ref:`Wang <pair-Wang>`, where a polyhedron does not have sharp corners
and edges, but is rounded at its vertices and edges by spheres
centered on each vertex with a specified diameter. The edges if the
centered on each vertex with a specified diameter. The edges of the
polyhedron are defined between pairs of adjacent vertices. Its faces
are defined by a loop of edges. The sphere diameter for each polygon
is specified in the data file read by the :doc:`read data <read_data>`

342
src/REPLICA/fix_pimd_langevin.cpp
View File

@ -50,7 +50,7 @@ using namespace FixConst;
using MathConst::MY_PI;
using MathConst::THIRD;
enum { NMPIMD };
enum { PIMD, NMPIMD };
enum { PHYSICAL, NORMAL };
enum { BAOAB, OBABO };
enum { ISO, ANISO, TRICLINIC };
@ -121,6 +121,8 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[i], "method") == 0) {
if (strcmp(arg[i + 1], "nmpimd") == 0)
method = NMPIMD;
else if (strcmp(arg[i + 1], "pimd") == 0)
method = PIMD;
else
error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin");
} else if (strcmp(arg[i], "integrator") == 0) {
@ -159,7 +161,7 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
error->universe_all(FLERR, "Invalid fmass value for fix pimd/langevin");
} else if (strcmp(arg[i], "sp") == 0) {
sp = utils::numeric(FLERR, arg[i + 1], false, lmp);
if (sp < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd/nvt");
if (sp < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd/langevin");
} else if (strcmp(arg[i], "fmmode") == 0) {
if (strcmp(arg[i + 1], "physical") == 0)
fmmode = PHYSICAL;
@ -170,9 +172,14 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
"Unknown fictitious mass mode for fix pimd/langevin. Only physical "
"mass and normal mode mass are supported!");
} else if (strcmp(arg[i], "scale") == 0) {
if (method == PIMD)
error->universe_all(
FLERR,
"The scale parameter of the PILE_L thermostat is not supported for method pimd. Delete "
"scale parameter if you do want to use method pimd.");
pilescale = utils::numeric(FLERR, arg[i + 1], false, lmp);
if (pilescale < 0.0)
error->universe_all(FLERR, "Invalid pile scale value for fix pimd/langevin");
error->universe_all(FLERR, "Invalid PILE_L scale value for fix pimd/langevin");
} else if (strcmp(arg[i], "temp") == 0) {
temp = utils::numeric(FLERR, arg[i + 1], false, lmp);
if (temp < 0.0) error->universe_all(FLERR, "Invalid temp value for fix pimd/langevin");
@ -245,6 +252,15 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
error->universe_all(
FLERR, fmt::format("Must not use pressure coupling with {} ensemble", Ensembles[ensemble]));
if (method == PIMD && pstat_flag)
error->universe_all(FLERR,
"Pressure control has not been supported for method pimd yet. Please set "
"method to nmpimd.");
if (method == PIMD && fmmode == NORMAL)
error->universe_all(
FLERR, "Normal mode mass is not supported for method pimd. Please set method to nmpimd.");
/* Initiation */
global_freq = 1;
@ -482,6 +498,13 @@ void FixPIMDLangevin::setup(int vflag)
nmpimd_transform(bufsortedall, x, M_x2xp[universe->iworld]);
else if (cmode == MULTI_PROC)
nmpimd_transform(bufbeads, x, M_x2xp[universe->iworld]);
} else if (method == PIMD) {
inter_replica_comm(x);
spring_force();
} else {
error->universe_all(
FLERR,
"Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
}
collect_xc();
compute_spring_energy();
@ -533,11 +556,18 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/)
nmpimd_transform(bufsortedall, x, M_x2xp[universe->iworld]);
else if (cmode == MULTI_PROC)
nmpimd_transform(bufbeads, x, M_x2xp[universe->iworld]);
qc_step();
a_step();
qc_step();
a_step();
} else if (method == PIMD) {
q_step();
q_step();
} else {
error->universe_all(
FLERR,
"Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
}
qc_step();
a_step();
qc_step();
a_step();
} else if (integrator == BAOAB) {
if (pstat_flag) {
compute_totke();
@ -551,23 +581,42 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/)
nmpimd_transform(bufsortedall, x, M_x2xp[universe->iworld]);
else if (cmode == MULTI_PROC)
nmpimd_transform(bufbeads, x, M_x2xp[universe->iworld]);
qc_step();
a_step();
} else if (method == PIMD) {
q_step();
} else {
error->universe_all(
FLERR,
"Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
}
qc_step();
a_step();
if (tstat_flag) {
o_step();
if (removecomflag) remove_com_motion();
if (pstat_flag) press_o_step();
}
qc_step();
a_step();
if (method == NMPIMD) {
qc_step();
a_step();
} else if (method == PIMD) {
q_step();
} else {
error->universe_all(
FLERR,
"Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
}
} else {
error->universe_all(FLERR, "Unknown integrator parameter for fix pimd/langevin");
error->universe_all(FLERR,
"Unknown integrator parameter for fix pimd/langevin. Only obabo and baoab "
"integrators are supported!");
}
collect_xc();
compute_spring_energy();
compute_t_prim();
compute_p_prim();
if (method == NMPIMD) {
compute_spring_energy();
compute_t_prim();
compute_p_prim();
}
if (method == NMPIMD) {
inter_replica_comm(x);
@ -609,30 +658,45 @@ void FixPIMDLangevin::final_integrate()
void FixPIMDLangevin::post_force(int /*flag*/)
{
if (atom->nmax > maxunwrap) reallocate_x_unwrap();
if (atom->nmax > maxxc) reallocate_xc();
int nlocal = atom->nlocal;
double **x = atom->x;
double **f = atom->f;
imageint *image = atom->image;
tagint *tag = atom->tag;
for (int i = 0; i < nlocal; i++) {
x_unwrap[i][0] = x[i][0];
x_unwrap[i][1] = x[i][1];
x_unwrap[i][2] = x[i][2];
}
if (mapflag) {
for (int i = 0; i < nlocal; i++) { domain->unmap(x_unwrap[i], image[i]); }
}
for (int i = 0; i < nlocal; i++) {
xc[i][0] = xcall[3 * (tag[i] - 1) + 0];
xc[i][1] = xcall[3 * (tag[i] - 1) + 1];
xc[i][2] = xcall[3 * (tag[i] - 1) + 2];
if (method == NMPIMD) {
if (atom->nmax > maxunwrap) reallocate_x_unwrap();
if (atom->nmax > maxxc) reallocate_xc();
for (int i = 0; i < nlocal; i++) {
x_unwrap[i][0] = x[i][0];
x_unwrap[i][1] = x[i][1];
x_unwrap[i][2] = x[i][2];
}
if (mapflag) {
for (int i = 0; i < nlocal; i++) { domain->unmap(x_unwrap[i], image[i]); }
}
for (int i = 0; i < nlocal; i++) {
xc[i][0] = xcall[3 * (tag[i] - 1) + 0];
xc[i][1] = xcall[3 * (tag[i] - 1) + 1];
xc[i][2] = xcall[3 * (tag[i] - 1) + 2];
}
compute_vir();
compute_cvir();
compute_t_vir();
}
compute_vir();
compute_cvir();
compute_t_vir();
if (method == PIMD) {
if (mapflag) {
for (int i = 0; i < nlocal; i++) { domain->unmap(x[i], image[i]); }
}
inter_replica_comm(x);
spring_force();
compute_spring_energy();
if (mapflag) {
for (int i = 0; i < nlocal; i++) { domain->unmap_inv(x[i], image[i]); }
}
}
compute_pote();
if (method == NMPIMD) {
inter_replica_comm(f);
@ -641,6 +705,7 @@ void FixPIMDLangevin::post_force(int /*flag*/)
else if (cmode == MULTI_PROC)
nmpimd_transform(bufbeads, f, M_x2xp[universe->iworld]);
}
c_pe->addstep(update->ntimestep + 1);
c_press->addstep(update->ntimestep + 1);
}
@ -655,6 +720,8 @@ void FixPIMDLangevin::end_of_step()
if (pstat_flag) compute_totenthalpy();
}
/* ---------------------------------------------------------------------- */
void FixPIMDLangevin::collect_xc()
{
int nlocal = atom->nlocal;
@ -688,7 +755,9 @@ void FixPIMDLangevin::collect_xc()
void FixPIMDLangevin::b_step()
{
// used for both NMPIMD and PIMD
// For NMPIMD, force only includes the contribution of external potential.
// For PIMD, force includes the contributions of external potential and spring force.
int n = atom->nlocal;
int *type = atom->type;
double **v = atom->v;
@ -706,6 +775,8 @@ void FixPIMDLangevin::b_step()
void FixPIMDLangevin::qc_step()
{
// used for NMPIMD
// evolve the centroid mode
int nlocal = atom->nlocal;
double **x = atom->x;
double **v = atom->v;
@ -770,6 +841,8 @@ void FixPIMDLangevin::qc_step()
void FixPIMDLangevin::a_step()
{
// used for NMPIMD
// use analytical solution of harmonic oscillator to evolve the non-centroid modes
int n = atom->nlocal;
double **x = atom->x;
double **v = atom->v;
@ -801,6 +874,25 @@ void FixPIMDLangevin::a_step()
/* ---------------------------------------------------------------------- */
void FixPIMDLangevin::q_step()
{
// used for PIMD
// evolve all beads
int nlocal = atom->nlocal;
double **x = atom->x;
double **v = atom->v;
if (!pstat_flag) {
for (int i = 0; i < nlocal; i++) {
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixPIMDLangevin::baro_init()
{
vw[0] = vw[1] = vw[2] = vw[3] = vw[4] = vw[5] = 0.0;
@ -895,20 +987,22 @@ void FixPIMDLangevin::langevin_init()
_omega_k = new double[np];
Lan_c = new double[np];
Lan_s = new double[np];
if (fmmode == PHYSICAL) {
for (int i = 0; i < np; i++) {
_omega_k[i] = _omega_np * sqrt(lam[i]) / sqrt(fmass);
Lan_c[i] = cos(sqrt(lam[i]) * _omega_np_dt_half);
Lan_s[i] = sin(sqrt(lam[i]) * _omega_np_dt_half);
if (method == NMPIMD) {
if (fmmode == PHYSICAL) {
for (int i = 0; i < np; i++) {
_omega_k[i] = _omega_np * sqrt(lam[i]) / sqrt(fmass);
Lan_c[i] = cos(sqrt(lam[i]) * _omega_np_dt_half);
Lan_s[i] = sin(sqrt(lam[i]) * _omega_np_dt_half);
}
} else if (fmmode == NORMAL) {
for (int i = 0; i < np; i++) {
_omega_k[i] = _omega_np / sqrt(fmass);
Lan_c[i] = cos(_omega_np_dt_half);
Lan_s[i] = sin(_omega_np_dt_half);
}
} else {
error->universe_all(FLERR, "Unknown fmmode setting; only physical and normal are supported!");
}
} else if (fmmode == NORMAL) {
for (int i = 0; i < np; i++) {
_omega_k[i] = _omega_np / sqrt(fmass);
Lan_c[i] = cos(_omega_np_dt_half);
Lan_s[i] = sin(_omega_np_dt_half);
}
} else {
error->universe_all(FLERR, "Unknown fmmode setting; only physical and normal are supported!");
}
if (tau > 0)
@ -930,27 +1024,34 @@ void FixPIMDLangevin::langevin_init()
if (thermostat == PILE_L) {
std::string out = "\nInitializing PI Langevin equation thermostat...\n";
out += "Bead ID | omega | tau | c1 | c2\n";
tau_k = new double[np];
c1_k = new double[np];
c2_k = new double[np];
tau_k[0] = tau;
c1_k[0] = c1;
c2_k[0] = c2;
for (int i = 1; i < np; i++) {
tau_k[i] = 0.5 / pilescale / _omega_k[i];
if (integrator == OBABO)
c1_k[i] = exp(-0.5 * update->dt / tau_k[i]);
else if (integrator == BAOAB)
c1_k[i] = exp(-1.0 * update->dt / tau_k[i]);
else
error->universe_all(FLERR,
"Unknown integrator parameter for fix pimd/langevin. Only obabo and "
"baoab integrators are supported!");
c2_k[i] = sqrt(1.0 - c1_k[i] * c1_k[i]);
}
for (int i = 0; i < np; i++) {
out += fmt::format(" {:d} {:.8e} {:.8e} {:.8e} {:.8e}\n", i, _omega_k[i], tau_k[i],
c1_k[i], c2_k[i]);
if (method == NMPIMD) {
tau_k = new double[np];
c1_k = new double[np];
c2_k = new double[np];
tau_k[0] = tau;
c1_k[0] = c1;
c2_k[0] = c2;
for (int i = 1; i < np; i++) {
tau_k[i] = 0.5 / pilescale / _omega_k[i];
if (integrator == OBABO)
c1_k[i] = exp(-0.5 * update->dt / tau_k[i]);
else if (integrator == BAOAB)
c1_k[i] = exp(-1.0 * update->dt / tau_k[i]);
else
error->universe_all(FLERR,
"Unknown integrator parameter for fix pimd/langevin. Only obabo and "
"baoab integrators are supported!");
c2_k[i] = sqrt(1.0 - c1_k[i] * c1_k[i]);
}
for (int i = 0; i < np; i++) {
out += fmt::format(" {:d} {:.8e} {:.8e} {:.8e} {:.8e}\n", i, _omega_k[i], tau_k[i],
c1_k[i], c2_k[i]);
}
} else if (method == PIMD) {
for (int i = 0; i < np; i++) {
out += fmt::format(" {:d} {:.8e} {:.8e} {:.8e} {:.8e}\n", i, _omega_np / sqrt(fmass),
tau, c1, c2);
}
}
if (thermostat == PILE_L) out += "PILE_L thermostat successfully initialized!\n";
out += "\n";
@ -966,13 +1067,24 @@ void FixPIMDLangevin::o_step()
int *type = atom->type;
double beta_np = 1.0 / force->boltz / Lan_temp * inverse_np * force->mvv2e;
if (thermostat == PILE_L) {
for (int i = 0; i < nlocal; i++) {
atom->v[i][0] = c1_k[universe->iworld] * atom->v[i][0] +
c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian();
atom->v[i][1] = c1_k[universe->iworld] * atom->v[i][1] +
c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian();
atom->v[i][2] = c1_k[universe->iworld] * atom->v[i][2] +
c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian();
if (method == NMPIMD) {
for (int i = 0; i < nlocal; i++) {
atom->v[i][0] = c1_k[universe->iworld] * atom->v[i][0] +
c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian();
atom->v[i][1] = c1_k[universe->iworld] * atom->v[i][1] +
c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian();
atom->v[i][2] = c1_k[universe->iworld] * atom->v[i][2] +
c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian();
}
} else if (method == PIMD) {
for (int i = 0; i < nlocal; i++) {
atom->v[i][0] =
c1 * atom->v[i][0] + c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian();
atom->v[i][1] =
c1 * atom->v[i][1] + c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian();
atom->v[i][2] =
c1 * atom->v[i][2] + c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian();
}
}
}
}
@ -1059,6 +1171,50 @@ void FixPIMDLangevin::nmpimd_transform(double **src, double **des, double *vecto
}
}
/* ---------------------------------------------------------------------- */
void FixPIMDLangevin::spring_force()
{
spring_energy = 0.0;
double **x = atom->x;
double **f = atom->f;
double *_mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
tagint *tagtmp = atom->tag;
// printf("iworld = %d, x_last = %d, x_next = %d\n", universe->iworld, x_last, x_next);
int *mask = atom->mask;
// int idx_tmp = atom->map(1);
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
double delx1 = bufsortedall[x_last * nlocal + tagtmp[i] - 1][0] - x[i][0];
double dely1 = bufsortedall[x_last * nlocal + tagtmp[i] - 1][1] - x[i][1];
double delz1 = bufsortedall[x_last * nlocal + tagtmp[i] - 1][2] - x[i][2];
double delx2 = bufsortedall[x_next * nlocal + tagtmp[i] - 1][0] - x[i][0];
double dely2 = bufsortedall[x_next * nlocal + tagtmp[i] - 1][1] - x[i][1];
double delz2 = bufsortedall[x_next * nlocal + tagtmp[i] - 1][2] - x[i][2];
double ff = fbond * _mass[type[i]];
// double ff = 0;
double dx = delx1 + delx2;
double dy = dely1 + dely2;
double dz = delz1 + delz2;
f[i][0] += (dx) *ff;
f[i][1] += (dy) *ff;
f[i][2] += (dz) *ff;
spring_energy += 0.5 * ff * (delx2 * delx2 + dely2 * dely2 + delz2 * delz2);
}
}
}
/* ----------------------------------------------------------------------
Comm operations
------------------------------------------------------------------------- */
@ -1084,6 +1240,9 @@ void FixPIMDLangevin::comm_init()
planrecv[i] = universe->root_proc[irecv];
modeindex[i] = irecv;
}
x_next = (universe->iworld + 1 + universe->nworlds) % (universe->nworlds);
x_last = (universe->iworld - 1 + universe->nworlds) % (universe->nworlds);
}
/* ---------------------------------------------------------------------- */
@ -1312,21 +1471,32 @@ void FixPIMDLangevin::compute_totke()
void FixPIMDLangevin::compute_spring_energy()
{
spring_energy = 0.0;
total_spring_energy = se_bead = 0.0;
if (method == NMPIMD) {
spring_energy = 0.0;
total_spring_energy = se_bead = 0.0;
double **x = atom->x;
double *_mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
double **x = atom->x;
double *_mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
spring_energy += 0.5 * _mass[type[i]] * fbond * lam[universe->iworld] *
(x[i][0] * x[i][0] + x[i][1] * x[i][1] + x[i][2] * x[i][2]);
for (int i = 0; i < nlocal; i++) {
spring_energy += 0.5 * _mass[type[i]] * fbond * lam[universe->iworld] *
(x[i][0] * x[i][0] + x[i][1] * x[i][1] + x[i][2] * x[i][2]);
}
MPI_Allreduce(&spring_energy, &se_bead, 1, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(&se_bead, &total_spring_energy, 1, MPI_DOUBLE, MPI_SUM, universe->uworld);
total_spring_energy /= universe->procs_per_world[universe->iworld];
} else if (method == PIMD) {
total_spring_energy = se_bead = 0.0;
MPI_Allreduce(&spring_energy, &se_bead, 1, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(&se_bead, &total_spring_energy, 1, MPI_DOUBLE, MPI_SUM, universe->uworld);
total_spring_energy /= universe->procs_per_world[universe->iworld];
} else {
error->universe_all(
FLERR,
"Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
}
MPI_Allreduce(&spring_energy, &se_bead, 1, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(&se_bead, &total_spring_energy, 1, MPI_DOUBLE, MPI_SUM, universe->uworld);
total_spring_energy /= universe->procs_per_world[universe->iworld];
}
/* ---------------------------------------------------------------------- */

4
src/REPLICA/fix_pimd_langevin.h
View File

@ -77,6 +77,8 @@ class FixPIMDLangevin : public Fix {
int me, nprocs, ireplica, nreplica, nprocs_universe;
int ntotal, maxlocal;
int x_last, x_next;
int cmode;
int sizeplan;
int *plansend, *planrecv;
@ -93,6 +95,7 @@ class FixPIMDLangevin : public Fix {
void comm_init();
void inter_replica_comm(double **ptr);
void spring_force();
/* normal-mode operations */
@ -121,6 +124,7 @@ class FixPIMDLangevin : public Fix {
a_step(); // integrate for dt/2 according to A part (non-centroid mode, harmonic force between replicas)
void qc_step(); // integrate for dt/2 for the centroid mode (x <- x + v * dt/2)
void o_step(); // integrate for dt according to O part (O-U process, for thermostating)
void q_step(); // integrate for dt/2 for all the beads (x <- x + v * dt/2)
/* Bussi-Zykova-Parrinello barostat */