Skip to content

Commit

Permalink
Merge pull request #472 from fweik/pmass
Browse files Browse the repository at this point in the history 
Removed PMASS macro
  • Loading branch information
@fweik
fweik committed Nov 30, 2015
2 parents 8948567 + 1937470 commit 6699f97
Show file tree
Hide file tree
Showing 20 changed files with 121 additions and 126 deletions.
4 changes: 2 additions & 2 deletions src/core/comfixed.cpp
View file
Expand Up @@ -88,7 +88,7 @@ void calc_comfixed()
np = cell->n;
for(i = 0; i < np; i++) {
if(p[i].p.type==t0) {
type_mass += PMASS(p[i]);
type_mass += (p[i]).p.mass;
for(j = 0; j < 3; j++) {
fsum0[j] += p[i].f.f[j];
}
Expand All @@ -103,7 +103,7 @@ void calc_comfixed()
for(i = 0; i < np; i++) {
if(p[i].p.type==t0) {
for(j = 0; j < 3; j++) {
p[i].f.f[j] -= fsum0[j]/type_mass*PMASS(p[i]);
p[i].f.f[j] -= fsum0[j]/type_mass*(p[i]).p.mass;
}
}
}
Expand Down
12 changes: 6 additions & 6 deletions src/core/communication.cpp
View file
Expand Up @@ -556,9 +556,9 @@ void mpi_send_f(int pnode, int part, double F[3])

if (pnode == this_node) {
Particle *p = local_particles[part];
F[0] /= PMASS(*p);
F[1] /= PMASS(*p);
F[2] /= PMASS(*p);
F[0] /= (*p).p.mass;
F[1] /= (*p).p.mass;
F[2] /= (*p).p.mass;
memmove(p->f.f, F, 3*sizeof(double));
}
else
Expand All @@ -573,9 +573,9 @@ void mpi_send_f_slave(int pnode, int part)
Particle *p = local_particles[part];
MPI_Recv(p->f.f, 3, MPI_DOUBLE, 0, SOME_TAG,
comm_cart, MPI_STATUS_IGNORE);
p->f.f[0] /= PMASS(*p);
p->f.f[1] /= PMASS(*p);
p->f.f[2] /= PMASS(*p);
p->f.f[0] /= (*p).p.mass;
p->f.f[1] /= (*p).p.mass;
p->f.f[2] /= (*p).p.mass;
}

on_particle_change();
Expand Down
8 changes: 4 additions & 4 deletions src/core/dpd.cpp
View file
Expand Up @@ -201,11 +201,11 @@ void add_dpd_thermo_pair_force(Particle * p1, Particle * p2, double d[3], double
#endif

#ifdef DPD_MASS_RED
massf=2*PMASS(*p1)*PMASS(*p2)/(PMASS(*p1)+PMASS(*p2));
massf=2*(*p1).p.mass*(*p2).p.mass/((*p1).p.mass+(*p2).p.mass);
#endif

#ifdef DPD_MASS_LIN
massf=0.5*(PMASS(*p1)+PMASS(*p2));
massf=0.5*((*p1).p.mass+(*p2).p.mass);
#endif


Expand Down Expand Up @@ -408,11 +408,11 @@ void add_inter_dpd_pair_force(Particle *p1, Particle *p2, IA_parameters *ia_para
if( le_chatterjee_test_pair(p1, p2) ) return;
#endif
#ifdef DPD_MASS_RED
massf=2*PMASS(*p1)*PMASS(*p2)/(PMASS(*p1)+PMASS(*p2));
massf=2*(*p1).p.mass*(*p2).p.mass/((*p1).p.mass+(*p2).p.mass);
#endif

#ifdef DPD_MASS_LIN
massf=0.5*(PMASS(*p1)+PMASS(*p2));
massf=0.5*((*p1).p.mass+(*p2).p.mass);
#endif

P_times_dist_sqr[0][0]=dist2;
Expand Down
4 changes: 2 additions & 2 deletions src/core/energy_inline.hpp
View file
Expand Up @@ -511,11 +511,11 @@ inline void add_kinetic_energy(Particle *p1)
// ostringstream msg;
// msg << "SMALL TIME STEP";
// energy.data.e[0] += SQR(smaller_time_step/time_step) *
// (SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*PMASS(*p1);
// (SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*(*p1).p.mass;
// }
// else
// #endif
energy.data.e[0] += (SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*PMASS(*p1);
energy.data.e[0] += (SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*(*p1).p.mass;

#ifdef ROTATION
#ifdef ROTATION_PER_PARTICLE
Expand Down
8 changes: 4 additions & 4 deletions src/core/ghmc.cpp
View file
Expand Up @@ -186,7 +186,7 @@ double calc_local_temp()
part = local_cells.cell[c]->part;
for (i = 0; i < np; i++) {
for (j = 0; j < 3; j++) {
temp += PMASS(part[i])*SQR(part[i].m.v[j]/time_step);
temp += (part[i]).p.mass*SQR(part[i].m.v[j]/time_step);
#ifdef ROTATION
temp += SQR(part[i].m.omega[j]);
#endif
Expand Down Expand Up @@ -219,7 +219,7 @@ void calc_kinetic(double *ek_trans , double *ek_rot)
#endif

/* kinetic energy */
et += (SQR(part[i].m.v[0]) + SQR(part[i].m.v[1]) + SQR(part[i].m.v[2]))*PMASS(part[i]);
et += (SQR(part[i].m.v[0]) + SQR(part[i].m.v[1]) + SQR(part[i].m.v[2]))*(part[i]).p.mass;

/* rotational energy */
#ifdef ROTATION
Expand Down Expand Up @@ -347,7 +347,7 @@ void simple_momentum_update()
#endif

#ifdef MASS
sigmat = sqrt(temperature / PMASS(part[i]));
sigmat = sqrt(temperature / (part[i]).p.mass);
#endif
for (j = 0; j < 3; j++) {
part[i].m.v[j] = sigmat*gaussian_random()*time_step;
Expand Down Expand Up @@ -380,7 +380,7 @@ void partial_momentum_update()
part = local_cells.cell[c]->part;
for (i = 0; i < np; i++) {
#ifdef MASS
sigmat = sqrt(temperature / PMASS(part[i]));
sigmat = sqrt(temperature / (part[i]).p.mass);
#endif
for (j = 0; j < 3; j++) {
part[i].m.v[j] = cosp*(part[i].m.v[j])+sinp*(sigmat*gaussian_random()*time_step);
Expand Down
20 changes: 10 additions & 10 deletions src/core/integrate.cpp
View file
Expand Up @@ -628,9 +628,9 @@ void rescale_forces()
np = cell->n;
for(i = 0; i < np; i++) {
check_particle_force(&p[i]);
p[i].f.f[0] *= scale/PMASS(p[i]);
p[i].f.f[1] *= scale/PMASS(p[i]);
p[i].f.f[2] *= scale/PMASS(p[i]);
p[i].f.f[0] *= scale/(p[i]).p.mass;
p[i].f.f[1] *= scale/(p[i]).p.mass;
p[i].f.f[2] *= scale/(p[i]).p.mass;

ONEPART_TRACE(if(p[i].p.identity==check_id) fprintf(stderr,"%d: OPT: SCAL f = (%.3e,%.3e,%.3e) v_old = (%.3e,%.3e,%.3e)\n",this_node,p[i].f.f[0],p[i].f.f[1],p[i].f.f[2],p[i].m.v[0],p[i].m.v[1],p[i].m.v[2]));

Expand Down Expand Up @@ -668,9 +668,9 @@ void rescale_forces_propagate_vel()
for(i = 0; i < np; i++) {
check_particle_force(&p[i]);
/* Rescale forces: f_rescaled = 0.5*dt*dt * f_calculated * (1/mass) */
p[i].f.f[0] *= scale/PMASS(p[i]);
p[i].f.f[1] *= scale/PMASS(p[i]);
p[i].f.f[2] *= scale/PMASS(p[i]);
p[i].f.f[0] *= scale/(p[i]).p.mass;
p[i].f.f[1] *= scale/(p[i]).p.mass;
p[i].f.f[2] *= scale/(p[i]).p.mass;

ONEPART_TRACE(if(p[i].p.identity==check_id) fprintf(stderr,"%d: OPT: SCAL f = (%.3e,%.3e,%.3e) v_old = (%.3e,%.3e,%.3e)\n",this_node,p[i].f.f[0],p[i].f.f[1],p[i].f.f[2],p[i].m.v[0],p[i].m.v[1],p[i].m.v[2]));
#ifdef VIRTUAL_SITES
Expand All @@ -683,13 +683,13 @@ void rescale_forces_propagate_vel()
#endif
#ifdef NPT
if(integ_switch == INTEG_METHOD_NPT_ISO && ( nptiso.geometry & nptiso.nptgeom_dir[j] )) {
nptiso.p_vel[j] += SQR(p[i].m.v[j])*PMASS(p[i]);
nptiso.p_vel[j] += SQR(p[i].m.v[j])*(p[i]).p.mass;
#ifdef MULTI_TIMESTEP
if (smaller_time_step > 0. && current_time_step_is_small==1)
p[i].m.v[j] += p[i].f.f[j];
else
#endif
p[i].m.v[j] += p[i].f.f[j] + friction_therm0_nptiso(p[i].m.v[j])/PMASS(p[i]);
p[i].m.v[j] += p[i].f.f[j] + friction_therm0_nptiso(p[i].m.v[j])/(p[i]).p.mass;
}
else
#endif
Expand Down Expand Up @@ -908,8 +908,8 @@ void propagate_vel()
p[i].m.v[j] += p[i].f.f[j];
else
#endif
p[i].m.v[j] += p[i].f.f[j] + friction_therm0_nptiso(p[i].m.v[j])/PMASS(p[i]);
nptiso.p_vel[j] += SQR(p[i].m.v[j])*PMASS(p[i]);
p[i].m.v[j] += p[i].f.f[j] + friction_therm0_nptiso(p[i].m.v[j])/(p[i]).p.mass;
nptiso.p_vel[j] += SQR(p[i].m.v[j])*(p[i]).p.mass;
}
else
#endif
Expand Down
4 changes: 2 additions & 2 deletions src/core/integrate_sd.cpp
View file
Expand Up @@ -448,7 +448,7 @@ void propagate_pos_sd()
for (int d=0;d<3;d++){
p[i].r.p[d] = pos[3*j+d]+box_l[d]*rint(p[i].r.p[d]/box_l[d]);
p[i].m.v[d] = velocity[3*j+d];
//p[i].f.f[d] *= (0.5*time_step*time_step)/PMASS(*part);
//p[i].f.f[d] *= (0.5*time_step*time_step)/(*part).p.mass;
}
#else
for (int d=0;d<3;d++){
Expand All @@ -458,7 +458,7 @@ void propagate_pos_sd()
#endif
// somehow this does not effect anything, although it is called ...
for (int d=0;d<3;d++){
p[i].f.f[d] *= (0.5*time_step*time_step)/PMASS(*p);
p[i].f.f[d] *= (0.5*time_step*time_step)/(*p).p.mass;
}
for (int d=0;d<DIM;d++){
assert(!isnan(pos[DIM*i+d]));
Expand Down
6 changes: 3 additions & 3 deletions src/core/molforces.cpp
View file
Expand Up @@ -202,14 +202,14 @@ void calc_local_mol_info (IntList *local_trapped_mols)
#endif
}
if (fixed) {
topology[mol].mass += PMASS(p[i]);
topology[mol].mass += (p[i]).p.mass;
/* Unfold the particle */
unfold_position(p[i].r.p, p[i].m.v, p[i].l.i);

for ( j = 0 ; j < 3 ; j++ ) {
topology[mol].f[j] += p[i].f.f[j];
topology[mol].com[j] += p[i].r.p[j]*PMASS(p[i]);
topology[mol].v[j] += p[i].m.v[j]*PMASS(p[i]);
topology[mol].com[j] += p[i].r.p[j]*(p[i]).p.mass;
topology[mol].v[j] += p[i].m.v[j]*(p[i]).p.mass;
}
/* Fold the particle back */
fold_position(p[i].r.p,p[i].l.i);
Expand Down
7 changes: 7 additions & 0 deletions src/core/particle_data.hpp
View file
Expand Up @@ -83,6 +83,13 @@ typedef struct {
#ifdef MASS
/** particle mass */
double mass;
#else
/** set mass to 0 */
#if HAVE_CXX11
constexpr static double mass = 1.0;
#else
const static double mass = 1.0;
#endif
#endif

#ifdef SHANCHEN
Expand Down
4 changes: 2 additions & 2 deletions src/core/pressure.cpp
View file
Expand Up @@ -976,8 +976,8 @@ int local_stress_tensor_calc(DoubleList *TensorInBin, int bins[3], int periodic[
PTENSOR_TRACE(fprintf(stderr,"%d:Ideal gas component is {",this_node));
for(k=0;k<3;k++) {
for(l=0;l<3;l++) {
TensorInBin[bin].e[k*3 + l] += (p1->m.v[k])*(p1->m.v[l])*PMASS(*p1)/time_step/time_step;
PTENSOR_TRACE(fprintf(stderr,"%f ",(p1->m.v[k])*(p1->m.v[l])*PMASS(*p1)/time_step/time_step));
TensorInBin[bin].e[k*3 + l] += (p1->m.v[k])*(p1->m.v[l])*(*p1).p.mass/time_step/time_step;
PTENSOR_TRACE(fprintf(stderr,"%f ",(p1->m.v[k])*(p1->m.v[l])*(*p1).p.mass/time_step/time_step));
}
}

Expand Down
12 changes: 6 additions & 6 deletions src/core/pressure.hpp
View file
Expand Up @@ -500,27 +500,27 @@ inline void add_kinetic_virials(Particle *p1,int v_comp)
#ifdef MULTI_TIMESTEP
if (smaller_time_step > 0.) {
if(v_comp)
virials.data.e[0] += SQR(time_step/smaller_time_step)*(SQR(p1->m.v[0] - p1->f.f[0]) + SQR(p1->m.v[1] - p1->f.f[1]) + SQR(p1->m.v[2] - p1->f.f[2]))*PMASS(*p1);
virials.data.e[0] += SQR(time_step/smaller_time_step)*(SQR(p1->m.v[0] - p1->f.f[0]) + SQR(p1->m.v[1] - p1->f.f[1]) + SQR(p1->m.v[2] - p1->f.f[2]))*(*p1).p.mass;
else
virials.data.e[0] += SQR(time_step/smaller_time_step)*(SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*PMASS(*p1);
virials.data.e[0] += SQR(time_step/smaller_time_step)*(SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*(*p1).p.mass;
} else
#endif
{
if(v_comp)
virials.data.e[0] += (SQR(p1->m.v[0] - p1->f.f[0]) + SQR(p1->m.v[1] - p1->f.f[1]) + SQR(p1->m.v[2] - p1->f.f[2]))*PMASS(*p1);
virials.data.e[0] += (SQR(p1->m.v[0] - p1->f.f[0]) + SQR(p1->m.v[1] - p1->f.f[1]) + SQR(p1->m.v[2] - p1->f.f[2]))*(*p1).p.mass;
else
virials.data.e[0] += (SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*PMASS(*p1);
virials.data.e[0] += (SQR(p1->m.v[0]) + SQR(p1->m.v[1]) + SQR(p1->m.v[2]))*(*p1).p.mass;
}

/* ideal gas contribution (the rescaling of the velocities by '/=time_step' each will be done later) */
for(k=0;k<3;k++)
for(l=0;l<3;l++)
#ifdef MULTI_TIMESTEP
if (smaller_time_step > 0.)
p_tensor.data.e[k*3 + l] += SQR(time_step/smaller_time_step)*(p1->m.v[k])*(p1->m.v[l])*PMASS(*p1);
p_tensor.data.e[k*3 + l] += SQR(time_step/smaller_time_step)*(p1->m.v[k])*(p1->m.v[l])*(*p1).p.mass;
else
#endif
p_tensor.data.e[k*3 + l] += (p1->m.v[k])*(p1->m.v[l])*PMASS(*p1);
p_tensor.data.e[k*3 + l] += (p1->m.v[k])*(p1->m.v[l])*(*p1).p.mass;

}

Expand Down
12 changes: 6 additions & 6 deletions src/core/rattle.cpp
View file
Expand Up @@ -160,14 +160,14 @@ void compute_pos_corr_vec(int *repeat_)
r_ij_dot = scalar(r_ij_t, r_ij);
G = 0.50*(ia_params->p.rigid_bond.d2 - r_ij2 )/r_ij_dot;
#ifdef MASS
G /= (PMASS(*p1)+PMASS(*p2));
G /= ((*p1).p.mass+(*p2).p.mass);
#else
G /= 2;
#endif
for (j=0;j<3;j++) {
pos_corr = G*r_ij_t[j];
p1->f.f[j] += pos_corr*PMASS(*p2);
p2->f.f[j] -= pos_corr*PMASS(*p1);
p1->f.f[j] += pos_corr*(*p2).p.mass;
p2->f.f[j] -= pos_corr*(*p1).p.mass;
}
/*Increase the 'repeat' flag by one */
*repeat_ = *repeat_ + 1;
Expand Down Expand Up @@ -313,15 +313,15 @@ void compute_vel_corr_vec(int *repeat_)
{
K = scalar(v_ij, r_ij)/ia_params->p.rigid_bond.d2;
#ifdef MASS
K /= (PMASS(*p1) + PMASS(*p2));
K /= ((*p1).p.mass + (*p2).p.mass);
#else
K /= 2.0;
#endif
for (j=0;j<3;j++)
{
vel_corr = K*r_ij[j];
p1->f.f[j] -= vel_corr*PMASS(*p2);
p2->f.f[j] += vel_corr*PMASS(*p1);
p1->f.f[j] -= vel_corr*(*p2).p.mass;
p2->f.f[j] += vel_corr*(*p1).p.mass;
}
*repeat_ = *repeat_ + 1 ;
}
Expand Down
12 changes: 6 additions & 6 deletions src/core/statistics.cpp
View file
Expand Up @@ -295,9 +295,9 @@ void centermass(int type, double *com)
for (j=0; j<n_part; j++) {
if ((partCfg[j].p.type == type) || (type == -1)) {
for (i=0; i<3; i++) {
com[i] += partCfg[j].r.p[i]*PMASS(partCfg[j]);
com[i] += partCfg[j].r.p[i]*(partCfg[j]).p.mass;
}
M += PMASS(partCfg[j]);
M += (partCfg[j]).p.mass;
}
}

Expand Down Expand Up @@ -343,7 +343,7 @@ void angularmomentum(int type, double *com)
if (type == partCfg[j].p.type)
{
vector_product(partCfg[j].r.p,partCfg[j].m.v,tmp);
pre_factor=PMASS(partCfg[j]);
pre_factor=(partCfg[j]).p.mass;
for (i=0; i<3; i++) {
com[i] += tmp[i]*pre_factor;
}
Expand All @@ -367,7 +367,7 @@ void momentofinertiamatrix(int type, double *MofImatrix)
for (i=0; i<3; i++) {
p1[i] = partCfg[j].r.p[i] - com[i];
}
massi= PMASS(partCfg[j]);
massi= (partCfg[j]).p.mass;
MofImatrix[0] += massi * (p1[1] * p1[1] + p1[2] * p1[2]) ;
MofImatrix[4] += massi * (p1[0] * p1[0] + p1[2] * p1[2]);
MofImatrix[8] += massi * (p1[0] * p1[0] + p1[1] * p1[1]);
Expand Down Expand Up @@ -1491,9 +1491,9 @@ void centermass_conf(int k, int type_1, double *com)
{
for (i=0; i<3; i++)
{
com[i] += configs[k][3*j+i]*PMASS(partCfg[j]);
com[i] += configs[k][3*j+i]*(partCfg[j]).p.mass;
}
M += PMASS(partCfg[j]);
M += (partCfg[j]).p.mass;
}
}
for (i=0; i<3; i++)
Expand Down

0 comments on commit 6699f97

Please sign in to comment.