Commit 99e5dc71 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

add support for fix_modify virial yes to fix smd

parent ccb67d8d
Loading
Loading
Loading
Loading
+5 −0
Original line number Diff line number Diff line
@@ -101,6 +101,11 @@ See the "read_restart"_read_restart.html command for info on how to
re-specify a fix in an input script that reads a restart file, so that
the operation of the fix continues in an uninterrupted fashion.

The "fix_modify"_fix_modify.html {virial} option is supported by this
fix to add the contribution due to the added forces on atoms to the
system's virial as part of "thermodynamic output"_thermo_style.html.
The default is {virial no}

The "fix_modify"_fix_modify.html {respa} option is supported by
this fix. This allows to set at which level of the "r-RESPA"_run_style.html
integrator the fix is adding its forces. Default is the outermost level.
+29 −0
Original line number Diff line number Diff line
@@ -58,6 +58,7 @@ FixSMD::FixSMD(LAMMPS *lmp, int narg, char **arg) :
  extvector = 1;
  respa_level_support = 1;
  ilevel_respa = 0;
  virial_flag = 1;

  int argoffs=3;
  if (strcmp(arg[argoffs],"cvel") == 0) {
@@ -181,6 +182,11 @@ void FixSMD::setup(int vflag)

void FixSMD::post_force(int vflag)
{
  // energy and virial setup

  if (vflag) v_setup(vflag);
  else evflag = 0;

  if (styleflag & SMD_TETHER) smd_tether();
  else smd_couple();

@@ -238,12 +244,15 @@ void FixSMD::smd_tether()
  // apply restoring force to atoms in group
  // f = -k*(r-r0)*mass/masstotal

  double **x = atom->x;
  double **f = atom->f;
  imageint *image = atom->image;
  int *mask = atom->mask;
  int *type = atom->type;
  double *mass = atom->mass;
  double *rmass = atom->rmass;
  double massfrac;
  double unwrap[3],v[6];
  int nlocal = atom->nlocal;

  ftotal[0] = ftotal[1] = ftotal[2] = 0.0;
@@ -259,6 +268,16 @@ void FixSMD::smd_tether()
        ftotal[0] -= fx*massfrac;
        ftotal[1] -= fy*massfrac;
        ftotal[2] -= fz*massfrac;
        if (evflag) {
          domain->unmap(x[i],image[i],unwrap);
          v[0] = fx*massfrac*unwrap[0];
          v[1] = fy*massfrac*unwrap[1];
          v[2] = fz*massfrac*unwrap[2];
          v[3] = fx*massfrac*unwrap[1];
          v[4] = fx*massfrac*unwrap[2];
          v[5] = fy*massfrac*unwrap[2];
          v_tally(i, v);
        }
      }
  } else {
    for (int i = 0; i < nlocal; i++)
@@ -270,6 +289,16 @@ void FixSMD::smd_tether()
        ftotal[0] -= fx*massfrac;
        ftotal[1] -= fy*massfrac;
        ftotal[2] -= fz*massfrac;
        if (evflag) {
          domain->unmap(x[i],image[i],unwrap);
          v[0] = fx*massfrac*unwrap[0];
          v[1] = fy*massfrac*unwrap[1];
          v[2] = fz*massfrac*unwrap[2];
          v[3] = fx*massfrac*unwrap[1];
          v[4] = fx*massfrac*unwrap[2];
          v[5] = fy*massfrac*unwrap[2];
          v_tally(i, v);
        }
      }
  }
}