diff --git a/src/set.cpp b/src/set.cpp index 37a0e815c4..56fbf42e54 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -430,8 +430,12 @@ void Set::command(int narg, char **arg) } else if (strcmp(arg[iarg],"dpd/theta") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); - if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1); - else dvalue = force->numeric(FLERR,arg[iarg+1]); + if (strcmp(arg[iarg+1],"NULL") == 0) dvalue = -1.0; + else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1); + else { + dvalue = force->numeric(FLERR,arg[iarg+1]); + if (dvalue < 0.0) error->all(FLERR,"Illegal set command"); + } if (!atom->dpd_flag) error->all(FLERR,"Cannot set dpd/theta for this atom style"); set(DPDTHETA); @@ -632,7 +636,20 @@ void Set::set(int keyword) atom->rmass[i] = atom->vfrac[i] * dvalue; } else if (keyword == SMD_CONTACT_RADIUS) atom->contact_radius[i] = dvalue; - else if (keyword == DPDTHETA) atom->dpdTheta[i] = dvalue; + + else if (keyword == DPDTHETA) { + if (dvalue >= 0.0) atom->dpdTheta[i] = dvalue; + else { + double onemass; + if (atom->rmass) onemass = atom->rmass[i]; + else onemass = atom->mass[atom->type[i]]; + double vx = atom->v[i][0]; + double vy = atom->v[i][1]; + double vz = atom->v[i][2]; + double tfactor = force->mvv2e / (domain->dimension * force->boltz); + atom->dpdTheta[i] = tfactor * onemass * (vx*vx + vy*vy + vz*vz); + } + } // set shape of ellipsoidal particle