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

First commit, things not working though.

Browse Source
This commit is contained in:
Sam Cameron
2022-03-29 07:50:19 +01:00
parent c30ba70fab
commit 0bebe6c428
7 changed files with 137 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/BROWNIAN/fix_brownian.cpp
View File

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

73
src/BROWNIAN/fix_brownian_asphere.cpp
View File

@ -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 <int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D>
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D, int Tp_2Drot>
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];

3
src/BROWNIAN/fix_brownian_asphere.h
View File

@ -35,8 +35,9 @@ class FixBrownianAsphere : public FixBrownianBase {
class AtomVecEllipsoid *avec;
private:
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D>
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D, int Tp_2Drot>
void initial_integrate_templated();
double g4;
};
} // namespace LAMMPS_NS
#endif

26
src/BROWNIAN/fix_brownian_base.cpp
View File

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

14
src/BROWNIAN/fix_brownian_base.h
View File

@ -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.
*/

48
src/BROWNIAN/fix_brownian_sphere.cpp
View File

@ -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 <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D>
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D, int Tp_2Drot>
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

7
src/BROWNIAN/fix_brownian_sphere.h
View File

@ -32,7 +32,8 @@ class FixBrownianSphere : public FixBrownianBase {
void initial_integrate(int);
private:
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D> void initial_integrate_templated();
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D, int Tp_2Drot>
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.