error->all(FLERR,"Cannot use fix qbmsst without per-type mass defined");
// set compute ptrs
int itemp = modify->find_compute(id_temp);
int ipress = modify->find_compute(id_press);
int ipe = modify->find_compute(id_pe);
if (itemp < 0 || ipress < 0|| ipe < 0)
error->all(FLERR,"Could not find fix qbmsst compute ID");
if (modify->compute[itemp]->tempflag == 0)
error->all(FLERR,"Fix qbmsst compute ID does not compute temperature");
if (modify->compute[ipress]->pressflag == 0)
error->all(FLERR,"Fix qbmsst compute ID does not compute pressure");
if (modify->compute[ipe]->peflag == 0)
error->all(FLERR,"Fix qbmsst compute ID does not compute potential energy");
temperature = modify->compute[itemp];
pressure = modify->compute[ipress];
pe = modify->compute[ipe];
// initiate the counter l and \mu
counter_l=0;
counter_mu=0;
// initiate qtb temperature
if (!qtb_set) {
t_current = t_init; qtb_set=1;
}
old_eavg = e0;
//set up the h time step for updating the random force \delta{}h=\frac{\pi}{\Omega_{max}}
if (int(1.0/(2*f_max*dtv)) == 0) {
if (comm->me == 0) printf ("Warning: Either f_max is too high or the time step is too big, setting f_max to be 1/timestep!\n");
h_timestep=dtv;
alpha=1;
} else {
alpha=int(1.0/(2*f_max*dtv));
h_timestep=alpha*dtv;
}
if (comm->me == 0) printf ("The effective maximum frequncy is now %f inverse time unit with alpha value as %d!\n", 0.5/h_timestep, alpha);
//gfactor is the random force \sqrt{\frac{2\gamma{}m_{i}}{\alpha*\delta{}t}}, \sqrt{12} makes the random array variance equal to unit.
for (int i = 1; i <= atom->ntypes; i++) {
gfactor[i] = sqrt(2*fric_coef*atom->mass[i])*sqrt(force->mvv2e)*sqrt(12/h_timestep);//this still leaves a square energy term from the power spectrum H.
}
// generate random number array with zero mean and variance equal 1/12.
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
for (int m = 0; m < 2*N_f; m++) {
random_array_0[i][m] = random->uniform()-0.5;
random_array_1[i][m] = random->uniform()-0.5;
random_array_2[i][m] = random->uniform()-0.5;
}
}
// initiate dilations
dilation[0] = dilation[1] = dilation[2] = 1.0;
// initialize the time derivative of the volume.
omega[0] = omega[1] = omega[2] = 0.0;
// compute total mass
double mass = 0.0;
for (int i = 0; i < atom->nlocal; i++) mass += atom->mass[atom->type[i]];