From 0bebe6c428b703f4f5aef7799a10aea231e23ea0 Mon Sep 17 00:00:00 2001 From: Sam Cameron Date: Tue, 29 Mar 2022 07:50:19 +0100 Subject: [PATCH] First commit, things not working though. --- src/BROWNIAN/fix_brownian.cpp | 5 +- src/BROWNIAN/fix_brownian_asphere.cpp | 73 +++++++++++++++++++-------- src/BROWNIAN/fix_brownian_asphere.h | 3 +- src/BROWNIAN/fix_brownian_base.cpp | 26 ++++++++-- src/BROWNIAN/fix_brownian_base.h | 14 ++++- src/BROWNIAN/fix_brownian_sphere.cpp | 48 ++++++++++++++---- src/BROWNIAN/fix_brownian_sphere.h | 7 ++- 7 files changed, 137 insertions(+), 39 deletions(-) diff --git a/src/BROWNIAN/fix_brownian.cpp b/src/BROWNIAN/fix_brownian.cpp index 4994b304e1..35afc17c87 100644 --- a/src/BROWNIAN/fix_brownian.cpp +++ b/src/BROWNIAN/fix_brownian.cpp @@ -33,7 +33,8 @@ using namespace FixConst; FixBrownian::FixBrownian(LAMMPS *lmp, int narg, char **arg) : FixBrownianBase(lmp, narg, arg) { - if (dipole_flag || gamma_t_eigen_flag || gamma_r_eigen_flag || gamma_r_flag) { + if (dipole_flag || gamma_t_eigen_flag || gamma_r_eigen_flag || gamma_r_flag + || rot_temp_flag || planar_rot_flag) { error->all(FLERR, "Illegal fix brownian command."); } if (!gamma_t_flag) { error->all(FLERR, "Illegal fix brownian command."); } @@ -45,7 +46,7 @@ void FixBrownian::init() { FixBrownianBase::init(); g1 /= gamma_t; - g2 *= sqrt(gamma_t); + g2 *= sqrt(temp/gamma_t); } /* ---------------------------------------------------------------------- */ diff --git a/src/BROWNIAN/fix_brownian_asphere.cpp b/src/BROWNIAN/fix_brownian_asphere.cpp index c2904aea64..7c65ecb8e9 100644 --- a/src/BROWNIAN/fix_brownian_asphere.cpp +++ b/src/BROWNIAN/fix_brownian_asphere.cpp @@ -85,6 +85,10 @@ void FixBrownianAsphere::init() } FixBrownianBase::init(); + + g4 = g2*sqrt(rot_temp); + g2 *= sqrt(temp); + } /* ---------------------------------------------------------------------- */ @@ -94,37 +98,55 @@ void FixBrownianAsphere::initial_integrate(int /*vflag */) if (domain->dimension == 2) { if (dipole_flag) { if (!noise_flag) { - initial_integrate_templated<0, 0, 1, 1>(); + initial_integrate_templated<0, 0, 1, 1, 0>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0, 1, 1, 1>(); + initial_integrate_templated<0, 1, 1, 1, 0>(); } else { - initial_integrate_templated<1, 0, 1, 1>(); + initial_integrate_templated<1, 0, 1, 1, 0>(); } } else { if (!noise_flag) { - initial_integrate_templated<0, 0, 0, 1>(); + initial_integrate_templated<0, 0, 0, 1, 0>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0, 1, 0, 1>(); + initial_integrate_templated<0, 1, 0, 1, 0>(); } else { - initial_integrate_templated<1, 0, 0, 1>(); + initial_integrate_templated<1, 0, 0, 1, 0>(); + } + } + } else if (planar_rot_flag) { + if (dipole_flag) { + if (!noise_flag) { + initial_integrate_templated<0, 0, 1, 0, 1>(); + } else if (gaussian_noise_flag) { + initial_integrate_templated<0, 1, 1, 0, 1>(); + } else { + initial_integrate_templated<1, 0, 1, 0, 1>(); + } + } else { + if (!noise_flag) { + initial_integrate_templated<0, 0, 0, 0, 1>(); + } else if (gaussian_noise_flag) { + initial_integrate_templated<0, 1, 0, 0, 1>(); + } else { + initial_integrate_templated<1, 0, 0, 0, 1>(); } } } else { if (dipole_flag) { if (!noise_flag) { - initial_integrate_templated<0, 0, 1, 0>(); + initial_integrate_templated<0, 0, 1, 0, 0>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0, 1, 1, 0>(); + initial_integrate_templated<0, 1, 1, 0, 0>(); } else { - initial_integrate_templated<1, 0, 1, 0>(); + initial_integrate_templated<1, 0, 1, 0, 0>(); } } else { if (!noise_flag) { - initial_integrate_templated<0, 0, 0, 0>(); + initial_integrate_templated<0, 0, 0, 0, 0>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0, 1, 0, 0>(); + initial_integrate_templated<0, 1, 0, 0, 0>(); } else { - initial_integrate_templated<1, 0, 0, 0>(); + initial_integrate_templated<1, 0, 0, 0, 0>(); } } } @@ -133,7 +155,7 @@ void FixBrownianAsphere::initial_integrate(int /*vflag */) /* ---------------------------------------------------------------------- */ -template +template void FixBrownianAsphere::initial_integrate_templated() { double **x = atom->x; @@ -172,21 +194,30 @@ void FixBrownianAsphere::initial_integrate_templated() if (Tp_2D) { tmp[0] = tmp[1] = 0.0; if (Tp_UNIFORM) { - tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g2; + tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g4; } else if (Tp_GAUSS) { - tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g2; + tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g4; + } else { + tmp[2] = g1 * tmp[2] * gamma_r_inv[2]; + } + } else if (Tp_2Drot) { + tmp[0] = tmp[1] = 0.0; + if (Tp_UNIFORM) { + tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g4; + } else if (Tp_GAUSS) { + tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g4; } else { tmp[2] = g1 * tmp[2] * gamma_r_inv[2]; } } else { if (Tp_UNIFORM) { - tmp[0] = g1 * tmp[0] * gamma_r_inv[0] + gamma_r_invsqrt[0] * (rng->uniform() - 0.5) * g2; - tmp[1] = g1 * tmp[1] * gamma_r_inv[1] + gamma_r_invsqrt[1] * (rng->uniform() - 0.5) * g2; - tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g2; + tmp[0] = g1 * tmp[0] * gamma_r_inv[0] + gamma_r_invsqrt[0] * (rng->uniform() - 0.5) * g4; + tmp[1] = g1 * tmp[1] * gamma_r_inv[1] + gamma_r_invsqrt[1] * (rng->uniform() - 0.5) * g4; + tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g4; } else if (Tp_GAUSS) { - tmp[0] = g1 * tmp[0] * gamma_r_inv[0] + gamma_r_invsqrt[0] * rng->gaussian() * g2; - tmp[1] = g1 * tmp[1] * gamma_r_inv[1] + gamma_r_invsqrt[1] * rng->gaussian() * g2; - tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g2; + tmp[0] = g1 * tmp[0] * gamma_r_inv[0] + gamma_r_invsqrt[0] * rng->gaussian() * g4; + tmp[1] = g1 * tmp[1] * gamma_r_inv[1] + gamma_r_invsqrt[1] * rng->gaussian() * g4; + tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g4; } else { tmp[0] = g1 * tmp[0] * gamma_r_inv[0]; tmp[1] = g1 * tmp[1] * gamma_r_inv[1]; diff --git a/src/BROWNIAN/fix_brownian_asphere.h b/src/BROWNIAN/fix_brownian_asphere.h index b40f797dd3..1a4fd54dfa 100644 --- a/src/BROWNIAN/fix_brownian_asphere.h +++ b/src/BROWNIAN/fix_brownian_asphere.h @@ -35,8 +35,9 @@ class FixBrownianAsphere : public FixBrownianBase { class AtomVecEllipsoid *avec; private: - template + template void initial_integrate_templated(); + double g4; }; } // namespace LAMMPS_NS #endif diff --git a/src/BROWNIAN/fix_brownian_base.cpp b/src/BROWNIAN/fix_brownian_base.cpp index 369d36bf70..6798b1f534 100644 --- a/src/BROWNIAN/fix_brownian_base.cpp +++ b/src/BROWNIAN/fix_brownian_base.cpp @@ -43,6 +43,8 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, n gamma_t_flag = gamma_r_flag = 0; gamma_t_eigen_flag = gamma_r_eigen_flag = 0; dipole_flag = 0; + rot_temp_flag = 0; + planar_rot_flag = 0; g2 = 0.0; if (narg < 5) error->all(FLERR, "Illegal fix brownian command."); @@ -154,11 +156,28 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, n if (gamma_r <= 0) error->all(FLERR, "Fix brownian gamma_r must be > 0."); iarg = iarg + 2; + } else if (strcmp(arg[iarg], "rotation_temp") == 0) { + if (narg == iarg + 1) { error->all(FLERR, "Illegal fix brownian command."); } + + rot_temp_flag = 1; + rot_temp = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (rot_temp <= 0) error->all(FLERR, "Fix brownian rotation_temp must be > 0."); + iarg = iarg + 2; + + } else if (strcmp(arg[iarg], "planar_rotation") == 0) { + + planar_rot_flag = 1; + if (domain->dimension == 2) + error->all(FLERR, "Do not explicitly set planar_rotation for 2D simulation"); + iarg = iarg + 1; + } else { error->all(FLERR, "Illegal fix brownian command."); } } - + if (!rot_temp_flag) { + rot_temp = temp; + } // initialize Marsaglia RNG with processor-unique seed rng = new RanMars(lmp, seed + comm->me); } @@ -196,14 +215,13 @@ void FixBrownianBase::init() { dt = update->dt; sqrtdt = sqrt(dt); - g1 = force->ftm2v; if (noise_flag == 0) { g2 = 0.0; } else if (gaussian_noise_flag == 1) { - g2 = sqrt(2 * force->boltz * temp / dt / force->mvv2e); + g2 = sqrt(2 * force->boltz / dt / force->mvv2e); } else { - g2 = sqrt(24 * force->boltz * temp / dt / force->mvv2e); + g2 = sqrt(24 * force->boltz / dt / force->mvv2e); } } diff --git a/src/BROWNIAN/fix_brownian_base.h b/src/BROWNIAN/fix_brownian_base.h index c13fc71a84..b1bb2b22e0 100644 --- a/src/BROWNIAN/fix_brownian_base.h +++ b/src/BROWNIAN/fix_brownian_base.h @@ -33,7 +33,9 @@ class FixBrownianBase : public Fix { int gamma_r_flag; // 0/1 if isotropic rotational damping is unset/set int gamma_t_eigen_flag; // 0/1 if anisotropic translational damping is unset/set int gamma_r_eigen_flag; // 0/1 if anisotropic rotational damping is unset/set - + int rot_temp_flag; // 0/1 if rotational temperature is unset/set + int planar_rot_flag; // 0/1 if rotation is constrained to 2D (xy) plane + double gamma_t, gamma_r; // translational and rotational (isotropic) damping params double *gamma_t_inv; // anisotropic damping parameter eigenvalues double *gamma_r_inv; @@ -47,6 +49,7 @@ class FixBrownianBase : public Fix { int gaussian_noise_flag; // 0/1 for uniform/gaussian noise double temp; // temperature + double rot_temp; // temperature double g1, g2; // prefactors in time stepping class RanMars *rng; @@ -85,4 +88,13 @@ E: Fix brownian gamma_r must be > 0. Self-explanatory. +E: Fix brownian rotation_temp must be > 0. + +Self-explanatory. + +E: Do not explicitly set planar_rotation for 2D simulation + +Self-explanatory. + + */ diff --git a/src/BROWNIAN/fix_brownian_sphere.cpp b/src/BROWNIAN/fix_brownian_sphere.cpp index c6664b18ef..292c0358c7 100644 --- a/src/BROWNIAN/fix_brownian_sphere.cpp +++ b/src/BROWNIAN/fix_brownian_sphere.cpp @@ -41,6 +41,9 @@ FixBrownianSphere::FixBrownianSphere(LAMMPS *lmp, int narg, char **arg) : if (!gamma_t_flag || !gamma_r_flag) { error->all(FLERR, "Illegal fix brownian command."); } if (!atom->mu_flag) error->all(FLERR, "Fix brownian/sphere requires atom attribute mu"); + if (!atom->sphere_flag) error->all(FLERR, "Fix brownian/sphere requires atom style sphere"); + + } /* ---------------------------------------------------------------------- */ @@ -50,9 +53,9 @@ void FixBrownianSphere::init() FixBrownianBase::init(); g3 = g1 / gamma_r; - g4 = g2 / sqrt(gamma_r); + g4 = g2 * sqrt(rot_temp/gamma_r); g1 /= gamma_t; - g2 /= sqrt(gamma_t); + g2 *= sqrt(temp/gamma_t); } /* ---------------------------------------------------------------------- */ @@ -61,19 +64,27 @@ void FixBrownianSphere::initial_integrate(int /*vflag */) { if (domain->dimension == 2) { if (!noise_flag) { - initial_integrate_templated<0, 0, 1>(); + initial_integrate_templated<0, 0, 1, 0>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0, 1, 1>(); + initial_integrate_templated<0, 1, 1, 0>(); } else { - initial_integrate_templated<1, 0, 1>(); + initial_integrate_templated<1, 0, 1, 0>(); + } + } else if (planar_rot_flag) { + if (!noise_flag) { + initial_integrate_templated<0, 0, 0, 1>(); + } else if (gaussian_noise_flag) { + initial_integrate_templated<0, 1, 0, 1>(); + } else { + initial_integrate_templated<1, 0, 0, 1>(); } } else { if (!noise_flag) { - initial_integrate_templated<0, 0, 0>(); + initial_integrate_templated<0, 0, 0,0>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0, 1, 0>(); + initial_integrate_templated<0, 1, 0,0>(); } else { - initial_integrate_templated<1, 0, 0>(); + initial_integrate_templated<1, 0, 0,0>(); } } return; @@ -81,7 +92,7 @@ void FixBrownianSphere::initial_integrate(int /*vflag */) /* ---------------------------------------------------------------------- */ -template +template void FixBrownianSphere::initial_integrate_templated() { double **x = atom->x; @@ -116,6 +127,24 @@ void FixBrownianSphere::initial_integrate_templated() dy = dt * g1 * f[i][1]; wz = 0; } + } else if (Tp_2Drot) { + wx = wy = 0; + if (Tp_UNIFORM) { + dx = dt * (g1 * f[i][0] + g2 * (rng->uniform() - 0.5)); + dy = dt * (g1 * f[i][1] + g2 * (rng->uniform() - 0.5)); + dz = dt * (g1 * f[i][2] + g2 * (rng->uniform() - 0.5)); + wz = (rng->uniform() - 0.5) * g4; + } else if (Tp_GAUSS) { + dx = dt * (g1 * f[i][0] + g2 * rng->gaussian()); + dy = dt * (g1 * f[i][1] + g2 * rng->gaussian()); + dz = dt * (g1 * f[i][2] + g2 * rng->gaussian()); + wz = rng->gaussian() * g4; + } else { + dx = dt * g1 * f[i][0]; + dy = dt * g1 * f[i][1]; + dz = dt * g1 * f[i][2]; + wz = 0; + } } else { if (Tp_UNIFORM) { dx = dt * (g1 * f[i][0] + g2 * (rng->uniform() - 0.5)); @@ -152,6 +181,7 @@ void FixBrownianSphere::initial_integrate_templated() wy += g3 * torque[i][1]; wz += g3 * torque[i][2]; + // store length of dipole as we need to convert it to a unit vector and // then back again diff --git a/src/BROWNIAN/fix_brownian_sphere.h b/src/BROWNIAN/fix_brownian_sphere.h index 73befb07f5..ae3b3b4b64 100644 --- a/src/BROWNIAN/fix_brownian_sphere.h +++ b/src/BROWNIAN/fix_brownian_sphere.h @@ -32,7 +32,8 @@ class FixBrownianSphere : public FixBrownianBase { void initial_integrate(int); private: - template void initial_integrate_templated(); + template + void initial_integrate_templated(); double g3, g4; }; } // namespace LAMMPS_NS @@ -53,6 +54,10 @@ E: Compute brownian/sphere requires atom attribute mu Self-explanatory. +E: Compute brownian/sphere requires atom style sphere + +Self-explanatory. + E: Fix brownian/sphere translational viscous drag coefficient must be > 0. Self-explanatory.