From 7550656ce7c424b9dd3263cd73239ec9412f6874 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 31 Dec 2015 11:57:15 -0500 Subject: [PATCH] add update dipole flag to nose-hoover integrator for extended particles --- doc/fix_nh.txt | 13 ++++++++++--- src/fix_nh.cpp | 13 +++++++++++++ src/fix_nh.h | 1 + src/fix_nh_sphere.cpp | 36 ++++++++++++++++++++++++++++++++++++ src/fix_nh_sphere.h | 1 + 5 files changed, 61 insertions(+), 3 deletions(-) diff --git a/doc/fix_nh.txt b/doc/fix_nh.txt index b288aa2e92..079c420ac3 100644 --- a/doc/fix_nh.txt +++ b/doc/fix_nh.txt @@ -26,7 +26,7 @@ fix ID group-ID style_name keyword value ... :pre ID, group-ID are documented in "fix"_fix.html command :ulb,l style_name = {nvt} or {npt} or {nph} :l one or more keyword/value pairs may be appended :l -keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {yz} or {xz} or {couple} or {tchain} or {pchain} or {mtk} or {tloop} or {ploop} or {nreset} or {drag} or {dilate} or {scalexy} or {scaleyz} or {scalexz} or {flip} or {fixedpoint} +keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {yz} or {xz} or {couple} or {tchain} or {pchain} or {mtk} or {tloop} or {ploop} or {nreset} or {drag} or {dilate} or {scalexy} or {scaleyz} or {scalexz} or {flip} or {fixedpoint} or {update} {temp} values = Tstart Tstop Tdamp Tstart,Tstop = external temperature at start/end of run Tdamp = temperature damping parameter (time units) @@ -56,8 +56,9 @@ keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {y {scalexz} value = {yes} or {no} = scale xz with lz {flip} value = {yes} or {no} = allow or disallow box flips when it becomes highly skewed {fixedpoint} values = x y z - x,y,z = perform barostat dilation/contraction around this point (distance units) :pre - + x,y,z = perform barostat dilation/contraction around this point (distance units) + {update} value = {dipole} update dipole orientation (only for sphere variants) :pre + :ule [Examples:] @@ -327,6 +328,12 @@ far. In all cases, the particle trajectories are unaffected by the chosen value, except for a time-dependent constant translation of positions. +If the {update} keyword is used with the {dipole} value, then the +orientation of the dipole moment of each particle is also updated +during the time integration. This option should be used for models +where a dipole moment is assigned to extended particles via use of +the "atom_style hybrid sphere dipole"_atom_style.html command. + :line NOTE: Using a barostat coupled to tilt dimensions {xy}, {xz}, {yz} can diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index 5f7bec6580..89bd562e8e 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -78,6 +78,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) omega_mass_flag = 0; etap_mass_flag = 0; flipflag = 1; + dipole_flag = 0; // turn on tilt factor scaling, whenever applicable @@ -321,6 +322,11 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; + } else if (strcmp(arg[iarg],"update") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); + if (strcmp(arg[iarg+1],"dipole") == 0) dipole_flag = 1; + else error->all(FLERR,"Illegal fix nvt/npt/nph command"); + iarg += 2; } else if (strcmp(arg[iarg],"fixedpoint") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]); @@ -417,6 +423,13 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) p_period[0] != p_period[2])) error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings"); + if (dipole_flag) { + if (!atom->sphere_flag) + error->all(FLERR,"Using update dipole flag requires atom style sphere"); + if (!atom->mu_flag) + error->all(FLERR,"Using update dipole flag requires atom attribute mu"); + } + if ((tstat_flag && t_period <= 0.0) || (p_flag[0] && p_period[0] <= 0.0) || (p_flag[1] && p_period[1] <= 0.0) || diff --git a/src/fix_nh.h b/src/fix_nh.h index be645d0be8..d2201e1612 100644 --- a/src/fix_nh.h +++ b/src/fix_nh.h @@ -109,6 +109,7 @@ class FixNH : public Fix { int eta_mass_flag; // 1 if eta_mass updated, 0 if not. int omega_mass_flag; // 1 if omega_mass updated, 0 if not. int etap_mass_flag; // 1 if etap_mass updated, 0 if not. + int dipole_flag; // 1 if dipole is updated, 0 if not. int scaleyz; // 1 if yz scaled with lz int scalexz; // 1 if xz scaled with lz diff --git a/src/fix_nh_sphere.cpp b/src/fix_nh_sphere.cpp index 9db1e34e39..f978b67660 100644 --- a/src/fix_nh_sphere.cpp +++ b/src/fix_nh_sphere.cpp @@ -91,6 +91,42 @@ void FixNHSphere::nve_v() } } +/* ---------------------------------------------------------------------- + perform full-step update of position with dipole orientation, if requested +-----------------------------------------------------------------------*/ + +void FixNHSphere::nve_x() +{ + // standard nve_x position update + + FixNH::nve_x(); + + // update mu for dipoles + // d_mu/dt = omega cross mu + // renormalize mu to dipole length + + if (dipole_flag) { + double msq,scale,g[3]; + double **mu = atom->mu; + double **omega = atom->omega; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (mu[i][3] > 0.0) { + g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]); + g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]); + g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]); + msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2]; + scale = mu[i][3]/sqrt(msq); + mu[i][0] = g[0]*scale; + mu[i][1] = g[1]*scale; + mu[i][2] = g[2]*scale; + } + } +} + /* ---------------------------------------------------------------------- perform half-step scaling of rotatonal velocities -----------------------------------------------------------------------*/ diff --git a/src/fix_nh_sphere.h b/src/fix_nh_sphere.h index a98356055d..0897f4a93b 100644 --- a/src/fix_nh_sphere.h +++ b/src/fix_nh_sphere.h @@ -26,6 +26,7 @@ class FixNHSphere : public FixNH { protected: void nve_v(); + void nve_x(); void nh_v_temp(); };