Hello, I have copied the update dipole code from the nve/sphere into the
void FixNHSphere::nve_v() at the end of the procedure (see attached), is
this enough, to make sure the calculation for dipoles is correct?
There is another point about "fix_nh", that irritated me, the temperture
calculation:
T = m * Sum(v^2) / dof, dof = degrees of freedom
In the case for dipole I would have 3N-3 translational and 2N rotational
dof (lammps uses 6N-3 by default).
Therefore I used T = 2 * N * (E_tot - E_pair) / dof, where N = number of
particles, with dof = 5N - 3. I also checked, if different temperature
definitions taking only into account E_rot or E_trans (see Allen &
Tildesley) result in the same temperature (after a time average), which is
the case in the "fix NVE/sphere update dipole" case.
T_trans = 2 * N * (E_tot - E_pair - E_rot) / (3N-3)
T_rot = 2 * N * (E_rot) / (2N)
Unfortunately in the NVT (Nose-Hover) T_trans was higher than T and T_rot
lower than T (significant in both cases) - what is the reason for this?
Does the temperature rescaling somehow affect this or does my update of
the dipoles fail and therefore the whole calculation includes errors?
By the way, a happy new year!
Best
G. Rosenthal
void FixNHSphere::nve_v()
{
// standard nve_v velocity update
FixNH::nve_v();
double **omega = atom->omega;
double **torque = atom->torque;
double *radius = atom->radius;
double *rmass = atom->rmass;
double *mass = atom->mass;
double **shape = atom->shape;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
// set timestep here since dt may have changed or come via rRESPA
double dtfrotate = dtf / INERTIA;
double dtirotate;
int itype;
// update omega for all particles
// d_omega/dt = torque / inertia
// 4 cases depending on radius vs shape and rmass vs mass
if (radius) {
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
omega[i][0] += dtirotate*torque[i][0];
omega[i][1] += dtirotate*torque[i][1];
omega[i][2] += dtirotate*torque[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtirotate = dtfrotate / (radius[i]*radius[i]*mass[type[i]]);
omega[i][0] += dtirotate*torque[i][0];
omega[i][1] += dtirotate*torque[i][1];
omega[i][2] += dtirotate*torque[i][2];
}
}
}
} else {
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
itype = type[i];
dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]);
omega[i][0] += dtirotate*torque[i][0];
omega[i][1] += dtirotate*torque[i][1];
omega[i][2] += dtirotate*torque[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
itype = type[i];
dtirotate = dtfrotate /
(shape[itype][0]*shape[itype][0]*mass[itype]);
omega[i][0] += dtirotate*torque[i][0];
omega[i][1] += dtirotate*torque[i][1];
omega[i][2] += dtirotate*torque[i][2];
}
}
}
}
// Begin modified
double g[3];
// update mu for dipoles
// d_mu/dt = omega cross mu
// renormalize mu to dipole length
double **mu = atom->mu;
double *dipole = atom->dipole;
double msq, scale;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (dipole[type[i]] > 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 = dipole[type[i]]/sqrt(msq);
mu[i][0] = g[0]*scale;
mu[i][1] = g[1]*scale;
mu[i][2] = g[2]*scale;
}
}
}
// End modified
}