tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

moved body momentum randomization into fix rigid small

Browse Source
This commit is contained in:
draneyj
2025-01-18 14:29:34 -05:00
parent 2e57a2ee23
commit d34a8a4807
4 changed files with 68 additions and 123 deletions
Download Patch File Download Diff File Expand all files Collapse all files

119
src/MC/fix_hmc.cpp
View File

@ -574,124 +574,7 @@ void FixHMC::random_velocities()
void FixHMC::rigid_body_random_velocities()
{
RigidSmallBody *body = fix_rigid->body;
int nlocal = fix_rigid->nlocal_body;
int ntotal = nlocal + fix_rigid->nghost_body;
int *mask = atom->mask;
double stdev, bmass, wbody[3], mbody[3];
double total_mass = 0;
RigidSmallBody *b;
double vcm[] = {0.0, 0.0, 0.0};
for (int ibody = 0; ibody < nlocal; ibody++) {
b = &body[ibody];
if (mask[b->ilocal] & groupbit) {
bmass = b->mass;
stdev = sqrt(KT / bmass);
total_mass += bmass;
for (int j = 0; j < 3; j++) {
b->vcm[j] = stdev * random->gaussian();
vcm[j] += b->vcm[j] * bmass;
if (b->inertia[j] > 0.0)
wbody[j] = sqrt(KT / b->inertia[j]) * random->gaussian();
else
wbody[j] = 0.0;
}
}
MathExtra::matvec(b->ex_space, b->ey_space, b->ez_space, wbody, b->omega);
}
if (mom_flag) {
for (int j = 0; j < 3; j++) vcm[j] /= total_mass;
for (int ibody = 0; ibody < nlocal; ibody++) {
if (mask[b->ilocal] & groupbit) {
b = &body[ibody];
for (int j = 0; j < 3; j++) b->vcm[j] -= vcm[j];
}
}
}
// Forward communicate vcm and omega to ghost bodies:
comm_flag = VCM_OMEGA;
comm->forward_comm(this, 6);
// Compute angular momenta of rigid bodies:
for (int ibody = 0; ibody < ntotal; ibody++) {
b = &body[ibody];
MathExtra::omega_to_angmom(b->omega, b->ex_space, b->ey_space, b->ez_space, b->inertia,
b->angmom);
MathExtra::transpose_matvec(b->ex_space, b->ey_space, b->ez_space, b->angmom, mbody);
MathExtra::quatvec(b->quat, mbody, b->conjqm);
b->conjqm[0] *= 2.0;
b->conjqm[1] *= 2.0;
b->conjqm[2] *= 2.0;
b->conjqm[3] *= 2.0;
}
// Compute velocities of individual atoms:
fix_rigid->set_v();
}
/* ----------------------------------------------------------------------
Pack rigid body info for forward communication
------------------------------------------------------------------------- */
int FixHMC::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
{
int *bodyown = fix_rigid->bodyown;
RigidSmallBody *body = fix_rigid->body;
int i, m, ibody;
RigidSmallBody *b;
m = 0;
for (i = 0; i < n; i++) {
ibody = bodyown[list[i]];
if (ibody >= 0) {
b = &body[ibody];
if (comm_flag == VCM_OMEGA) {
memcpy(&buf[m], b->vcm, three);
memcpy(&buf[m + 3], b->omega, three);
m += 6;
} else if (comm_flag == XCM) {
memcpy(&buf[m], b->xcm, three);
m += 3;
} else if (comm_flag == ROTATION) {
memcpy(&buf[m], b->ex_space, three);
memcpy(&buf[m + 3], b->ey_space, three);
memcpy(&buf[m + 6], b->ez_space, three);
memcpy(&buf[m + 9], b->quat, four);
m += 12;
}
}
}
return m;
}
/* ----------------------------------------------------------------------
Unpack rigid body info from forward communication
------------------------------------------------------------------------- */
void FixHMC::unpack_forward_comm(int n, int first, double *buf)
{
int *bodyown = fix_rigid->bodyown;
RigidSmallBody *body = fix_rigid->body;
int i, m, last;
RigidSmallBody *b;
m = 0;
last = first + n;
for (i = first; i < last; i++)
if (bodyown[i] >= 0) {
b = &body[bodyown[i]];
if (comm_flag == VCM_OMEGA) {
memcpy(b->vcm, &buf[m], three);
memcpy(b->omega, &buf[m + 3], three);
m += 6;
}
}
fix_rigid->resample_momenta(groupbit, mom_flag, random, KT);
}
/* ----------------------------------------------------------------------

3
src/MC/fix_hmc.h
View File

@ -35,9 +35,6 @@ class FixHMC : public Fix {
void end_of_step() override;
double compute_scalar() override;
double compute_vector(int) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
double memory_usage() override;
private:

66
src/RIGID/fix_rigid_small.cpp
View File

@ -33,6 +33,7 @@
#include "molecule.h"
#include "neighbor.h"
#include "random_mars.h"
#include "random_park.h"
#include "respa.h"
#include "rigid_const.h"
#include "tokenizer.h"
@ -2708,6 +2709,69 @@ void FixRigidSmall::write_restart_file(const char *file)
if (comm->me == 0) fclose(fp);
}
/* ----------------------------------------------------------------------
randomize velocities with respect to a given temperature
------------------------------------------------------------------------- */
void FixRigidSmall::resample_momenta(int groupbit, int mom_flag, class RanPark * random, double KT) {
int nlocal = nlocal_body;
int ntotal = nghost_body;
int *mask = atom->mask;
double stdev, bmass, wbody[3], mbody[3];
double total_mass = 0;
RigidSmallBody *b;
double vcm[] = {0.0, 0.0, 0.0};
for (int ibody = 0; ibody < nlocal; ibody++) {
b = &body[ibody];
if (mask[b->ilocal] & groupbit) {
bmass = b->mass;
stdev = sqrt(KT / bmass);
total_mass += bmass;
for (int j = 0; j < 3; j++) {
b->vcm[j] = stdev * random->gaussian();
vcm[j] += b->vcm[j] * bmass;
if (b->inertia[j] > 0.0)
wbody[j] = sqrt(KT / b->inertia[j]) * random->gaussian();
else
wbody[j] = 0.0;
}
}
MathExtra::matvec(b->ex_space, b->ey_space, b->ez_space, wbody, b->omega);
}
if (mom_flag) {
for (int j = 0; j < 3; j++) vcm[j] /= total_mass;
for (int ibody = 0; ibody < nlocal; ibody++) {
if (mask[b->ilocal] & groupbit) {
b = &body[ibody];
for (int j = 0; j < 3; j++) b->vcm[j] -= vcm[j];
}
}
}
// Forward communicate vcm and omega to ghost bodies:
commflag = FINAL;
comm->forward_comm(this, 10);
// Compute angular momenta of rigid bodies:
for (int ibody = 0; ibody < ntotal; ibody++) {
b = &body[ibody];
MathExtra::omega_to_angmom(b->omega, b->ex_space, b->ey_space, b->ez_space, b->inertia,
b->angmom);
MathExtra::transpose_matvec(b->ex_space, b->ey_space, b->ez_space, b->angmom, mbody);
MathExtra::quatvec(b->quat, mbody, b->conjqm);
b->conjqm[0] *= 2.0;
b->conjqm[1] *= 2.0;
b->conjqm[2] *= 2.0;
b->conjqm[3] *= 2.0;
}
// Compute velocities of individual atoms:
set_v();
}
/* ----------------------------------------------------------------------
allocate local atom-based arrays
------------------------------------------------------------------------- */
@ -3716,4 +3780,4 @@ void FixRigidSmall::check(int flag)
}
}
}
*/
*/

3
src/RIGID/fix_rigid_small.h
View File

@ -63,6 +63,7 @@ class FixRigidSmall : public Fix {
void initial_integrate_respa(int, int, int) override;
void final_integrate_respa(int, int) override;
void write_restart_file(const char *) override;
void resample_momenta(int, int, class RanPark *, double);
void grow_arrays(int) override;
void copy_arrays(int, int, int) override;
@ -219,4 +220,4 @@ class FixRigidSmall : public Fix {
} // namespace LAMMPS_NS
#endif
#endif
#endif