forked from lijiext/lammps
Merge branch 'collected-small-changes' into collected-post-stable-patches
This commit is contained in:
commit
fd18660463
|
@ -56,13 +56,13 @@ void FixEnforce2DKokkos<DeviceType>::post_force(int vflag)
|
|||
v = atomKK->k_v.view<DeviceType>();
|
||||
f = atomKK->k_f.view<DeviceType>();
|
||||
|
||||
if( atomKK->omega_flag )
|
||||
if (atomKK->omega_flag)
|
||||
omega = atomKK->k_omega.view<DeviceType>();
|
||||
|
||||
if( atomKK->angmom_flag )
|
||||
if (atomKK->angmom_flag)
|
||||
angmom = atomKK->k_angmom.view<DeviceType>();
|
||||
|
||||
if( atomKK->torque_flag )
|
||||
if (atomKK->torque_flag)
|
||||
torque = atomKK->k_torque.view<DeviceType>();
|
||||
|
||||
|
||||
|
@ -72,9 +72,9 @@ void FixEnforce2DKokkos<DeviceType>::post_force(int vflag)
|
|||
if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst;
|
||||
|
||||
int flag_mask = 0;
|
||||
if( atomKK->omega_flag ) flag_mask |= 1;
|
||||
if( atomKK->angmom_flag ) flag_mask |= 2;
|
||||
if( atomKK->torque_flag ) flag_mask |= 4;
|
||||
if (atomKK->omega_flag) flag_mask |= 1;
|
||||
if (atomKK->angmom_flag) flag_mask |= 2;
|
||||
if (atomKK->torque_flag) flag_mask |= 4;
|
||||
|
||||
copymode = 1;
|
||||
switch( flag_mask ){
|
||||
|
|
|
@ -1115,7 +1115,7 @@ int FixQEqReaxKokkos<DeviceType>::pack_forward_comm(int n, int *list, double *bu
|
|||
{
|
||||
int m;
|
||||
|
||||
if( pack_flag == 1)
|
||||
if (pack_flag == 1)
|
||||
for(m = 0; m < n; m++) buf[m] = h_d[list[m]];
|
||||
else if( pack_flag == 2 )
|
||||
for(m = 0; m < n; m++) buf[m] = h_s[list[m]];
|
||||
|
@ -1134,7 +1134,7 @@ void FixQEqReaxKokkos<DeviceType>::unpack_forward_comm(int n, int first, double
|
|||
{
|
||||
int i, m;
|
||||
|
||||
if( pack_flag == 1)
|
||||
if (pack_flag == 1)
|
||||
for(m = 0, i = first; m < n; m++, i++) h_d[i] = buf[m];
|
||||
else if( pack_flag == 2)
|
||||
for(m = 0, i = first; m < n; m++, i++) h_s[i] = buf[m];
|
||||
|
|
|
@ -554,7 +554,7 @@ void PairReaxCKokkos<DeviceType>::Deallocate_Lookup_Tables()
|
|||
|
||||
for( i = 0; i <= ntypes; ++i ) {
|
||||
for( j = i; j <= ntypes; ++j )
|
||||
if( LR[i][j].n ) {
|
||||
if (LR[i][j].n) {
|
||||
sfree( LR[i][j].y, "LR[i,j].y" );
|
||||
sfree( LR[i][j].H, "LR[i,j].H" );
|
||||
sfree( LR[i][j].vdW, "LR[i,j].vdW" );
|
||||
|
@ -1317,7 +1317,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeTabulatedLJCoulomb<N
|
|||
|
||||
/* Cubic Spline Interpolation */
|
||||
int r = (int)(rij * t.inv_dx);
|
||||
if( r == 0 ) ++r;
|
||||
if (r == 0) ++r;
|
||||
const F_FLOAT base = (double)(r+1) * t.dx;
|
||||
const F_FLOAT dif = rij - base;
|
||||
|
||||
|
@ -1534,7 +1534,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxBuildListsFull, const int &
|
|||
// hbond list
|
||||
if (i < nlocal && cut_hbsq > 0.0 && (ihb == 1 || ihb == 2) && rsq <= cut_hbsq) {
|
||||
jhb = paramssing(jtype).p_hbond;
|
||||
if( ihb == 1 && jhb == 2) {
|
||||
if (ihb == 1 && jhb == 2) {
|
||||
const int jj_index = hb_index - hb_first_i;
|
||||
if (jj_index >= maxhb) {
|
||||
d_resize_hb() = 1;
|
||||
|
@ -1702,7 +1702,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxBuildListsHalf<NEIGHFLAG>,
|
|||
// hbond list
|
||||
if (i < nlocal && cut_hbsq > 0.0 && (ihb == 1 || ihb == 2) && rsq <= cut_hbsq) {
|
||||
jhb = paramssing(jtype).p_hbond;
|
||||
if( ihb == 1 && jhb == 2) {
|
||||
if (ihb == 1 && jhb == 2) {
|
||||
if (NEIGHFLAG == HALF) {
|
||||
j_index = hb_first_i + d_hb_num[i];
|
||||
d_hb_num[i]++;
|
||||
|
@ -1925,7 +1925,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxBuildListsHalf_LessAtomics<
|
|||
// hbond list
|
||||
if (i < nlocal && cut_hbsq > 0.0 && (ihb == 1 || ihb == 2) && rsq <= cut_hbsq) {
|
||||
jhb = paramssing(jtype).p_hbond;
|
||||
if( ihb == 1 && jhb == 2) {
|
||||
if (ihb == 1 && jhb == 2) {
|
||||
if (NEIGHFLAG == HALF) {
|
||||
j_index = hb_first_i + d_hb_num[i];
|
||||
d_hb_num[i]++;
|
||||
|
@ -2288,7 +2288,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxBondOrder3, const int &ii)
|
|||
const F_FLOAT Clp = 2.0 * gp[15] * explp1 * (2.0 + vlpex);
|
||||
d_dDelta_lp[i] = Clp;
|
||||
|
||||
if( paramssing(itype).mass > 21.0 ) {
|
||||
if (paramssing(itype).mass > 21.0) {
|
||||
nlp_temp = 0.5 * (paramssing(itype).valency_e - paramssing(itype).valency);
|
||||
d_Delta_lp_temp[i] = paramssing(itype).nlp_opt - nlp_temp;
|
||||
} else {
|
||||
|
@ -2578,13 +2578,13 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeAngular<NEIGHFLAG,EV
|
|||
SBO = SBOp + (1.0 - prod_SBO) * (-d_Delta_boc[i] - p_val8 * vlpadj);
|
||||
dSBO1 = -8.0 * prod_SBO * ( d_Delta_boc[i] + p_val8 * vlpadj );
|
||||
|
||||
if( SBO <= 0.0 ) {
|
||||
if (SBO <= 0.0) {
|
||||
SBO2 = 0.0;
|
||||
CSBO2 = 0.0;
|
||||
} else if( SBO > 0.0 && SBO <= 1.0 ) {
|
||||
} else if (SBO > 0.0 && SBO <= 1.0) {
|
||||
SBO2 = pow( SBO, p_val9 );
|
||||
CSBO2 = p_val9 * pow( SBO, p_val9 - 1.0 );
|
||||
} else if( SBO > 1.0 && SBO < 2.0 ) {
|
||||
} else if (SBO > 1.0 && SBO < 2.0) {
|
||||
SBO2 = 2.0 - pow( 2.0-SBO, p_val9 );
|
||||
CSBO2 = p_val9 * pow( 2.0 - SBO, p_val9 - 1.0 );
|
||||
} else {
|
||||
|
@ -2640,8 +2640,8 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeAngular<NEIGHFLAG,EV
|
|||
// theta and derivatives
|
||||
|
||||
cos_theta = (delij[0]*delik[0]+delij[1]*delik[1]+delij[2]*delik[2])/(rij*rik);
|
||||
if( cos_theta > 1.0 ) cos_theta = 1.0;
|
||||
if( cos_theta < -1.0 ) cos_theta = -1.0;
|
||||
if (cos_theta > 1.0) cos_theta = 1.0;
|
||||
if (cos_theta < -1.0) cos_theta = -1.0;
|
||||
theta = acos(cos_theta);
|
||||
|
||||
const F_FLOAT inv_dists = 1.0 / (rij * rik);
|
||||
|
@ -2682,7 +2682,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeAngular<NEIGHFLAG,EV
|
|||
theta_0 = theta_0*constPI/180.0;
|
||||
|
||||
expval2theta = exp( -p_val2 * (theta_0-theta)*(theta_0-theta) );
|
||||
if( p_val1 >= 0 )
|
||||
if (p_val1 >= 0)
|
||||
expval12theta = p_val1 * (1.0 - expval2theta);
|
||||
else // To avoid linear Me-H-Me angles (6/6/06)
|
||||
expval12theta = p_val1 * -expval2theta;
|
||||
|
@ -2910,8 +2910,8 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeTorsion<NEIGHFLAG,EV
|
|||
const F_FLOAT rik = sqrt(rsqik);
|
||||
|
||||
cos_ijk = (delij[0]*delik[0]+delij[1]*delik[1]+delij[2]*delik[2])/(rij*rik);
|
||||
if( cos_ijk > 1.0 ) cos_ijk = 1.0;
|
||||
if( cos_ijk < -1.0 ) cos_ijk = -1.0;
|
||||
if (cos_ijk > 1.0) cos_ijk = 1.0;
|
||||
if (cos_ijk < -1.0) cos_ijk = -1.0;
|
||||
theta_ijk = acos(cos_ijk);
|
||||
|
||||
// dcos_ijk
|
||||
|
@ -2925,7 +2925,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeTorsion<NEIGHFLAG,EV
|
|||
}
|
||||
|
||||
sin_ijk = sin( theta_ijk );
|
||||
if( sin_ijk >= 0 && sin_ijk <= 1e-10 )
|
||||
if (sin_ijk >= 0 && sin_ijk <= 1e-10)
|
||||
tan_ijk_i = cos_ijk / 1e-10;
|
||||
else if( sin_ijk <= 0 && sin_ijk >= -1e-10 )
|
||||
tan_ijk_i = -cos_ijk / 1e-10;
|
||||
|
@ -2952,8 +2952,8 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeTorsion<NEIGHFLAG,EV
|
|||
BOA_jl = bo_jl - thb_cut;
|
||||
|
||||
cos_jil = -(delij[0]*deljl[0]+delij[1]*deljl[1]+delij[2]*deljl[2])/(rij*rjl);
|
||||
if( cos_jil > 1.0 ) cos_jil = 1.0;
|
||||
if( cos_jil < -1.0 ) cos_jil = -1.0;
|
||||
if (cos_jil > 1.0) cos_jil = 1.0;
|
||||
if (cos_jil < -1.0) cos_jil = -1.0;
|
||||
theta_jil = acos(cos_jil);
|
||||
|
||||
// dcos_jil
|
||||
|
@ -2968,7 +2968,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeTorsion<NEIGHFLAG,EV
|
|||
}
|
||||
|
||||
sin_jil = sin( theta_jil );
|
||||
if( sin_jil >= 0 && sin_jil <= 1e-10 )
|
||||
if (sin_jil >= 0 && sin_jil <= 1e-10)
|
||||
tan_jil_i = cos_jil / 1e-10;
|
||||
else if( sin_jil <= 0 && sin_jil >= -1e-10 )
|
||||
tan_jil_i = -cos_jil / 1e-10;
|
||||
|
@ -3008,21 +3008,21 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeTorsion<NEIGHFLAG,EV
|
|||
hnhe = rik * rjl * sin_ijk * cos_jil;
|
||||
|
||||
poem = 2.0 * rik * rjl * sin_ijk * sin_jil;
|
||||
if( poem < 1e-20 ) poem = 1e-20;
|
||||
if (poem < 1e-20) poem = 1e-20;
|
||||
|
||||
tel = SQR(rik) + SQR(rij) + SQR(rjl) - SQR(rlk) -
|
||||
2.0 * (rik * rij * cos_ijk - rik * rjl * cos_ijk * cos_jil + rij * rjl * cos_jil);
|
||||
|
||||
arg = tel / poem;
|
||||
if( arg > 1.0 ) arg = 1.0;
|
||||
if( arg < -1.0 ) arg = -1.0;
|
||||
if (arg > 1.0) arg = 1.0;
|
||||
if (arg < -1.0) arg = -1.0;
|
||||
|
||||
F_FLOAT sin_ijk_rnd = sin_ijk;
|
||||
F_FLOAT sin_jil_rnd = sin_jil;
|
||||
|
||||
if( sin_ijk >= 0 && sin_ijk <= 1e-10 ) sin_ijk_rnd = 1e-10;
|
||||
if (sin_ijk >= 0 && sin_ijk <= 1e-10) sin_ijk_rnd = 1e-10;
|
||||
else if( sin_ijk <= 0 && sin_ijk >= -1e-10 ) sin_ijk_rnd = -1e-10;
|
||||
if( sin_jil >= 0 && sin_jil <= 1e-10 ) sin_jil_rnd = 1e-10;
|
||||
if (sin_jil >= 0 && sin_jil <= 1e-10) sin_jil_rnd = 1e-10;
|
||||
else if( sin_jil <= 0 && sin_jil >= -1e-10 ) sin_jil_rnd = -1e-10;
|
||||
|
||||
// dcos_omega_di
|
||||
|
@ -3196,7 +3196,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeHydrogen<NEIGHFLAG,E
|
|||
|
||||
const int i = d_ilist[ii];
|
||||
const int itype = type(i);
|
||||
if( paramssing(itype).p_hbond != 1 ) return;
|
||||
if (paramssing(itype).p_hbond != 1) return;
|
||||
|
||||
const X_FLOAT xtmp = x(i,0);
|
||||
const X_FLOAT ytmp = x(i,1);
|
||||
|
@ -3216,7 +3216,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeHydrogen<NEIGHFLAG,E
|
|||
const int j_index = jj - j_start;
|
||||
const F_FLOAT bo_ij = d_BO(i,j_index);
|
||||
|
||||
if( paramssing(jtype).p_hbond == 2 && bo_ij >= HB_THRESHOLD ) {
|
||||
if (paramssing(jtype).p_hbond == 2 && bo_ij >= HB_THRESHOLD) {
|
||||
hblist[top] = jj;
|
||||
top ++;
|
||||
}
|
||||
|
@ -3256,8 +3256,8 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeHydrogen<NEIGHFLAG,E
|
|||
|
||||
// theta and derivatives
|
||||
cos_theta = (delij[0]*delik[0]+delij[1]*delik[1]+delij[2]*delik[2])/(rij*rik);
|
||||
if( cos_theta > 1.0 ) cos_theta = 1.0;
|
||||
if( cos_theta < -1.0 ) cos_theta = -1.0;
|
||||
if (cos_theta > 1.0) cos_theta = 1.0;
|
||||
if (cos_theta < -1.0) cos_theta = -1.0;
|
||||
theta = acos(cos_theta);
|
||||
|
||||
const F_FLOAT inv_dists = 1.0 / (rij * rik);
|
||||
|
@ -3475,8 +3475,8 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeBond1<NEIGHFLAG,EVFL
|
|||
// Stabilisation terminal triple bond
|
||||
F_FLOAT estriph = 0.0;
|
||||
|
||||
if( BO_i >= 1.00 ) {
|
||||
if( gp[37] == 2 || (imass == 12.0000 && jmass == 15.9990) ||
|
||||
if (BO_i >= 1.00) {
|
||||
if (gp[37] == 2 || (imass == 12.0000 && jmass == 15.9990) ||
|
||||
(jmass == 12.0000 && imass == 15.9990) ) {
|
||||
const F_FLOAT exphu = exp(-gp[7] * SQR(BO_i - 2.50) );
|
||||
const F_FLOAT exphua1 = exp(-gp[3] * (d_total_bo[i]-BO_i));
|
||||
|
|
|
@ -208,7 +208,7 @@ void PairComb3::coeff(int narg, char **arg)
|
|||
nelements = 0;
|
||||
for (i = 3; i < narg; i++) {
|
||||
if ((strcmp(arg[i],"C") == 0) && (cflag == 0)) {
|
||||
if( comm->me == 0 && screen) fprintf(screen,
|
||||
if (comm->me == 0 && screen) fprintf(screen,
|
||||
" PairComb3: Found C: reading additional library file\n");
|
||||
read_lib();
|
||||
cflag = 1;
|
||||
|
@ -922,7 +922,7 @@ void PairComb3::Short_neigh()
|
|||
|
||||
icontrol = params[iparam_ij].jelementgp;
|
||||
|
||||
if( icontrol == 1)
|
||||
if (icontrol == 1)
|
||||
xcctmp[i] += comb_fc(rr1,¶ms[iparam_ij]) * params[iparam_ij].pcross;
|
||||
if (icontrol == 2)
|
||||
xchtmp[i] += comb_fc(rr1,¶ms[iparam_ij]) * params[iparam_ij].pcross;
|
||||
|
@ -1250,7 +1250,7 @@ void PairComb3::compute(int eflag, int vflag)
|
|||
|
||||
// torsion: i-j-k-l: apply to all C-C bonds
|
||||
|
||||
if( params[iparam_ij].tor_flag != 0 ) {
|
||||
if (params[iparam_ij].tor_flag != 0) {
|
||||
srmu = vec3_dot(delrj,delrk)/(sqrt(rsq1*rsq2));
|
||||
srmu = sqrt(1.0-srmu*srmu);
|
||||
|
||||
|
@ -1382,7 +1382,7 @@ void PairComb3::compute(int eflag, int vflag)
|
|||
}
|
||||
|
||||
// torsion and radical: apply to all C-C bonds
|
||||
if( params[iparam_ijk].tor_flag != 0 && fabs(ptorr)>1.0e-8) {
|
||||
if (params[iparam_ijk].tor_flag != 0 && fabs(ptorr)>1.0e-8) {
|
||||
srmu = vec3_dot(delrj,delrk)/(sqrt(rsq1*rsq2));
|
||||
srmu = sqrt(1.0-srmu*srmu);
|
||||
|
||||
|
@ -2578,7 +2578,7 @@ void PairComb3::tables()
|
|||
rvdw[1][inty] = params[iparam_ij].vsig * 0.950;
|
||||
|
||||
// radius check: outer radius vs. sigma
|
||||
if( rvdw[0][inty] > rvdw[1][inty] )
|
||||
if (rvdw[0][inty] > rvdw[1][inty])
|
||||
error->all(FLERR,"Error in vdw spline: inner radius > outer radius");
|
||||
|
||||
rrc[0] = rvdw[1][inty];
|
||||
|
|
|
@ -122,7 +122,7 @@ void PairLCBOP::allocate()
|
|||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLCBOP::settings(int narg, char **/*arg*/) {
|
||||
if( narg != 0 ) error->all(FLERR,"Illegal pair_style command");
|
||||
if (narg != 0 ) error->all(FLERR,"Illegal pair_style command");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
@ -402,7 +402,7 @@ void PairLCBOP::FSR(int eflag, int /*vflag*/)
|
|||
r_sq = delx*delx + dely*dely + delz*delz;
|
||||
rijmag = sqrt(r_sq);
|
||||
f_c_ij = f_c( rijmag,r_1,r_2,&df_c_ij );
|
||||
if( f_c_ij <= TOL ) continue;
|
||||
if (f_c_ij <= TOL) continue;
|
||||
|
||||
VR = A*exp(-alpha*rijmag);
|
||||
dVRdi = -alpha*VR;
|
||||
|
@ -502,10 +502,10 @@ void PairLCBOP::FLR(int eflag, int /*vflag*/)
|
|||
df_c_ij = -df_c_ij;
|
||||
// derivative may be inherited from previous call, see f_c_LR definition
|
||||
f_c_ij *= f_c_LR( rijmag, r_1_LR, r_2_LR, &df_c_ij );
|
||||
if( f_c_ij <= TOL ) continue;
|
||||
if (f_c_ij <= TOL) continue;
|
||||
|
||||
V = dVdi = 0;
|
||||
if( rijmag<r_0 ) {
|
||||
if (rijmag<r_0) {
|
||||
double exp_part = exp( -lambda_1*(rijmag-r_0) );
|
||||
V = eps_1*( exp_part*exp_part - 2*exp_part) + v_1;
|
||||
dVdi = 2*eps_1*lambda_1*exp_part*( 1-exp_part );
|
||||
|
@ -553,7 +553,7 @@ void PairLCBOP::FNij( int i, int j, double factor, double **f, int vflag_atom )
|
|||
rik[1] = x[atomi][1]-x[atomk][1];
|
||||
rik[2] = x[atomi][2]-x[atomk][2];
|
||||
double riksq = (rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2]);
|
||||
if( riksq > r_1*r_1 ) { // && riksq < r_2*r_2, if second condition not fulfilled neighbor would not be in the list
|
||||
if (riksq > r_1*r_1) { // && riksq < r_2*r_2, if second condition not fulfilled neighbor would not be in the list
|
||||
double rikmag = sqrt(riksq);
|
||||
double df_c_ik;
|
||||
f_c( rikmag, r_1, r_2, &df_c_ik );
|
||||
|
@ -598,7 +598,7 @@ void PairLCBOP::FMij( int i, int j, double factor, double **f, int vflag_atom )
|
|||
double Fx = 1-f_c_LR(Nki, 2,3,&dF);
|
||||
dF = -dF;
|
||||
|
||||
if( df_c_ik > TOL ) {
|
||||
if (df_c_ik > TOL) {
|
||||
double factor2 = factor*df_c_ik*Fx;
|
||||
// F = factor2*(-grad rikmag)
|
||||
// grad_i rikmag = \vec{rik} /rikmag
|
||||
|
@ -613,7 +613,7 @@ void PairLCBOP::FMij( int i, int j, double factor, double **f, int vflag_atom )
|
|||
if (vflag_atom) v_tally2(atomi,atomk,fpair,rik);
|
||||
}
|
||||
|
||||
if( dF > TOL ) {
|
||||
if (dF > TOL) {
|
||||
double factor2 = factor*f_c_ik*dF;
|
||||
FNij( atomk, atomi, factor2, f, vflag_atom );
|
||||
}
|
||||
|
@ -676,12 +676,12 @@ double PairLCBOP::bondorder(int i, int j, double rij[3],
|
|||
double num_Nconj = ( Nij+1 )*( Nji+1 )*( Nij_el+Nji_el ) - 4*( Nij+Nji+2);
|
||||
double den_Nconj = Nij*( 3-Nij )*( Nji+1 ) + Nji*( 3-Nji )*( Nij+1 ) + eps;
|
||||
Nconj = num_Nconj / den_Nconj;
|
||||
if( Nconj <= 0 ) {
|
||||
if (Nconj <= 0) {
|
||||
Nconj = 0;
|
||||
dNconj_dNij = 0;
|
||||
dNconj_dNji = 0;
|
||||
dNconj_dNel = 0;
|
||||
} else if( Nconj >= 1 ) {
|
||||
} else if (Nconj >= 1) {
|
||||
Nconj = 1;
|
||||
dNconj_dNij = 0;
|
||||
dNconj_dNji = 0;
|
||||
|
@ -703,21 +703,21 @@ double PairLCBOP::bondorder(int i, int j, double rij[3],
|
|||
Fij_conj = F_conj( Nij, Nji, Nconj, &dF_dNij, &dF_dNji, &dF_dNconj );
|
||||
|
||||
/*forces for Nij*/
|
||||
if( 3-Nij > TOL ) {
|
||||
if (3-Nij > TOL) {
|
||||
double factor = -VA*0.5*( dF_dNij + dF_dNconj*( dNconj_dNij + dNconj_dNel*dNij_el_dNij ) );
|
||||
FNij( i, j, factor, f, vflag_atom );
|
||||
}
|
||||
/*forces for Nji*/
|
||||
if( 3-Nji > TOL ) {
|
||||
if (3-Nji > TOL) {
|
||||
double factor = -VA*0.5*( dF_dNji + dF_dNconj*( dNconj_dNji + dNconj_dNel*dNji_el_dNji ) );
|
||||
FNij( j, i, factor, f, vflag_atom );
|
||||
}
|
||||
/*forces for Mij*/
|
||||
if( 3-Mij > TOL ) {
|
||||
if (3-Mij > TOL) {
|
||||
double factor = -VA*0.5*( dF_dNconj*dNconj_dNel*dNij_el_dMij );
|
||||
FMij( i, j, factor, f, vflag_atom );
|
||||
}
|
||||
if( 3-Mji > TOL ) {
|
||||
if (3-Mji > TOL) {
|
||||
double factor = -VA*0.5*( dF_dNconj*dNconj_dNel*dNji_el_dMji );
|
||||
FMij( j, i, factor, f, vflag_atom );
|
||||
}
|
||||
|
@ -900,14 +900,14 @@ double PairLCBOP::gSpline( double x, double *dgdc ) {
|
|||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLCBOP::hSpline( double x, double *dhdx ) {
|
||||
if( x < -d ) {
|
||||
if (x < -d) {
|
||||
double z = kappa*( x+d );
|
||||
double y = pow(z, 10.0);
|
||||
double w = pow( 1+y, -0.1 );
|
||||
*dhdx = kappa*L*w/(1+y);
|
||||
return L*( 1 + z*w );
|
||||
}
|
||||
if( x > d ) {
|
||||
if (x > d) {
|
||||
*dhdx = R_1;
|
||||
return R_0 + R_1*( x-d );
|
||||
}
|
||||
|
@ -941,13 +941,13 @@ double PairLCBOP::F_conj( double N_ij, double N_ji, double N_conj_ij, double *dF
|
|||
double dF_0_dx = 0, dF_0_dy = 0;
|
||||
double dF_1_dx = 0, dF_1_dy = 0;
|
||||
double l, r;
|
||||
if( N_conj_ij < 1 ) {
|
||||
if (N_conj_ij < 1) {
|
||||
l = (1-y)* (1-x); r = ( f0.f_00 + x* x* f0.f_x_10 + y* y* f0.f_y_01 ); F_0 += l*r; dF_0_dx += -(1-y)*r +l*2*x* f0.f_x_10; dF_0_dy += -(1-x)*r +l*2*y* f0.f_y_01;
|
||||
l = (1-y)* x; r = ( f0.f_10 + (1-x)*(1-x)*f0.f_x_00 + y* y* f0.f_y_11 ); F_0 += l*r; dF_0_dx += (1-y)*r -l*2*(1-x)*f0.f_x_00; dF_0_dy += -x* r +l*2*y* f0.f_y_11;
|
||||
l = y* (1-x); r = ( f0.f_01 + x* x* f0.f_x_11 + (1-y)*(1-y)*f0.f_y_00 ); F_0 += l*r; dF_0_dx += -y* r +l*2*x* f0.f_x_11; dF_0_dy += (1-x)*r -l*2*(1-y)*f0.f_y_00;
|
||||
l = y* x; r = ( f0.f_11 + (1-x)*(1-x)*f0.f_x_01 + (1-y)*(1-y)*f0.f_y_10 ); F_0 += l*r; dF_0_dx += y* r -l*2*(1-x)*f0.f_x_01; dF_0_dy += x* r -l*2*(1-y)*f0.f_y_10;
|
||||
}
|
||||
if( N_conj_ij > 0 ) {
|
||||
if (N_conj_ij > 0) {
|
||||
l = (1-y)* (1-x); r = ( f0.f_00 + x* x* f1.f_x_10 + y* y* f1.f_y_01 ); F_1 += l*r; dF_1_dx += -(1-y)*r +l*2*x* f1.f_x_10; dF_1_dy += -(1-x)*r +l*2*y* f1.f_y_01;
|
||||
l = (1-y)* x; r = ( f1.f_10 + (1-x)*(1-x)*f1.f_x_00 + y* y* f1.f_y_11 ); F_1 += l*r; dF_1_dx += (1-y)*r -l*2*(1-x)*f1.f_x_00; dF_1_dy += -x* r +l*2*y* f1.f_y_11;
|
||||
l = y* (1-x); r = ( f1.f_01 + x* x* f1.f_x_11 + (1-y)*(1-y)*f1.f_y_00 ); F_1 += l*r; dF_1_dx += -y* r +l*2*x* f1.f_x_11; dF_1_dy += (1-x)*r -l*2*(1-y)*f1.f_y_00;
|
||||
|
|
|
@ -70,7 +70,7 @@ static xdr_uint32_t xdr_swapbytes(xdr_uint32_t x)
|
|||
static xdr_uint32_t xdr_htonl(xdr_uint32_t x)
|
||||
{
|
||||
short s=0x0F00;
|
||||
if( *((char *)&s)==(char)0x0F) {
|
||||
if (*((char *)&s)==(char)0x0F) {
|
||||
/* bigendian, do nothing */
|
||||
return x;
|
||||
} else {
|
||||
|
@ -82,7 +82,7 @@ static xdr_uint32_t xdr_htonl(xdr_uint32_t x)
|
|||
static xdr_uint32_t xdr_ntohl(xdr_uint32_t x)
|
||||
{
|
||||
short s=0x0F00;
|
||||
if( *((char *)&s)==(char)0x0F) {
|
||||
if (*((char *)&s)==(char)0x0F) {
|
||||
/* bigendian, do nothing */
|
||||
return x;
|
||||
} else {
|
||||
|
|
|
@ -204,7 +204,7 @@ double BondGromos::single(int type, double rsq, int /*i*/, int /*j*/,
|
|||
void *BondGromos::extract( char *str, int &dim )
|
||||
{
|
||||
dim = 1;
|
||||
if( strcmp(str,"kappa")==0) return (void*) k;
|
||||
if( strcmp(str,"r0")==0) return (void*) r0;
|
||||
if (strcmp(str,"kappa")==0) return (void*) k;
|
||||
if (strcmp(str,"r0")==0) return (void*) r0;
|
||||
return NULL;
|
||||
}
|
||||
|
|
|
@ -207,8 +207,8 @@ double BondHarmonic::single(int type, double rsq, int /*i*/, int /*j*/,
|
|||
void *BondHarmonic::extract( char *str, int &dim )
|
||||
{
|
||||
dim = 1;
|
||||
if( strcmp(str,"kappa")==0) return (void*) k;
|
||||
if( strcmp(str,"r0")==0) return (void*) r0;
|
||||
if (strcmp(str,"kappa")==0) return (void*) k;
|
||||
if (strcmp(str,"r0")==0) return (void*) r0;
|
||||
return NULL;
|
||||
}
|
||||
|
||||
|
|
|
@ -475,7 +475,7 @@ int FixQEq::pack_forward_comm(int n, int *list, double *buf,
|
|||
{
|
||||
int m;
|
||||
|
||||
if( pack_flag == 1)
|
||||
if (pack_flag == 1)
|
||||
for(m = 0; m < n; m++) buf[m] = d[list[m]];
|
||||
else if( pack_flag == 2 )
|
||||
for(m = 0; m < n; m++) buf[m] = s[list[m]];
|
||||
|
@ -493,7 +493,7 @@ void FixQEq::unpack_forward_comm(int n, int first, double *buf)
|
|||
{
|
||||
int i, m;
|
||||
|
||||
if( pack_flag == 1)
|
||||
if (pack_flag == 1)
|
||||
for(m = 0, i = first; m < n; m++, i++) d[i] = buf[m];
|
||||
else if( pack_flag == 2)
|
||||
for(m = 0, i = first; m < n; m++, i++) s[i] = buf[m];
|
||||
|
|
|
@ -103,7 +103,7 @@ void FixQEqDynamic::pre_force(int /*vflag*/)
|
|||
|
||||
if (update->ntimestep % nevery) return;
|
||||
|
||||
if( atom->nmax > nmax ) reallocate_storage();
|
||||
if (atom->nmax > nmax) reallocate_storage();
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
|
@ -251,7 +251,7 @@ int FixQEqDynamic::pack_forward_comm(int n, int *list, double *buf,
|
|||
{
|
||||
int m=0;
|
||||
|
||||
if( pack_flag == 1 )
|
||||
if (pack_flag == 1)
|
||||
for(m = 0; m < n; m++) buf[m] = atom->q[list[m]];
|
||||
else if( pack_flag == 2 )
|
||||
for(m = 0; m < n; m++) buf[m] = qf[list[m]];
|
||||
|
@ -265,7 +265,7 @@ void FixQEqDynamic::unpack_forward_comm(int n, int first, double *buf)
|
|||
{
|
||||
int i, m;
|
||||
|
||||
if( pack_flag == 1)
|
||||
if (pack_flag == 1)
|
||||
for(m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m];
|
||||
else if( pack_flag == 2)
|
||||
for(m = 0, i = first; m < n; m++, i++) qf[i] = buf[m];
|
||||
|
|
|
@ -122,7 +122,7 @@ void FixQEqFire::pre_force(int /*vflag*/)
|
|||
|
||||
if (ntimestep % nevery) return;
|
||||
|
||||
if( atom->nmax > nmax ) reallocate_storage();
|
||||
if (atom->nmax > nmax) reallocate_storage();
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
|
@ -315,7 +315,7 @@ int FixQEqFire::pack_forward_comm(int n, int *list, double *buf,
|
|||
{
|
||||
int m = 0;
|
||||
|
||||
if( pack_flag == 1 )
|
||||
if (pack_flag == 1)
|
||||
for(m = 0; m < n; m++) buf[m] = atom->q[list[m]];
|
||||
else if( pack_flag == 2 )
|
||||
for(m = 0; m < n; m++) buf[m] = qf[list[m]];
|
||||
|
@ -329,7 +329,7 @@ void FixQEqFire::unpack_forward_comm(int n, int first, double *buf)
|
|||
{
|
||||
int i, m;
|
||||
|
||||
if( pack_flag == 1)
|
||||
if (pack_flag == 1)
|
||||
for(m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m];
|
||||
else if( pack_flag == 2)
|
||||
for(m = 0, i = first; m < n; m++, i++) qf[i] = buf[m];
|
||||
|
|
|
@ -73,9 +73,9 @@ void FixQEqPoint::pre_force(int /*vflag*/)
|
|||
|
||||
nlocal = atom->nlocal;
|
||||
|
||||
if( atom->nmax > nmax ) reallocate_storage();
|
||||
if (atom->nmax > nmax) reallocate_storage();
|
||||
|
||||
if( nlocal > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE )
|
||||
if (nlocal > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE)
|
||||
reallocate_matrix();
|
||||
|
||||
init_matvec();
|
||||
|
|
|
@ -117,9 +117,9 @@ void FixQEqShielded::pre_force(int /*vflag*/)
|
|||
|
||||
nlocal = atom->nlocal;
|
||||
|
||||
if( atom->nmax > nmax ) reallocate_storage();
|
||||
if (atom->nmax > nmax) reallocate_storage();
|
||||
|
||||
if( nlocal > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE )
|
||||
if (nlocal > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE)
|
||||
reallocate_matrix();
|
||||
|
||||
init_matvec();
|
||||
|
|
|
@ -114,9 +114,9 @@ void FixQEqSlater::pre_force(int /*vflag*/)
|
|||
nlocal = atom->nlocal;
|
||||
nall = atom->nlocal + atom->nghost;
|
||||
|
||||
if( atom->nmax > nmax ) reallocate_storage();
|
||||
if (atom->nmax > nmax) reallocate_storage();
|
||||
|
||||
if( nlocal > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE )
|
||||
if (nlocal > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE)
|
||||
reallocate_matrix();
|
||||
|
||||
init_matvec();
|
||||
|
|
|
@ -147,12 +147,12 @@ int NBinSSA::coord2ssaAIR(const double *x)
|
|||
if(iz < 0){
|
||||
return -1;
|
||||
} else if(iz == 0){
|
||||
if( iy<0 ) return -1; // bottom left/middle/right
|
||||
if( (iy==0) && (ix<0) ) return -1; // left atoms
|
||||
if( (iy==0) && (ix==0) ) return 0; // Locally owned atoms
|
||||
if( (iy==0) && (ix>0) ) return 2; // Right atoms
|
||||
if( (iy>0) && (ix==0) ) return 1; // Top-middle atoms
|
||||
if( (iy>0) && (ix!=0) ) return 3; // Top-right and top-left atoms
|
||||
if (iy<0) return -1; // bottom left/middle/right
|
||||
if ((iy==0) && (ix<0) ) return -1; // left atoms
|
||||
if ((iy==0) && (ix==0)) return 0; // Locally owned atoms
|
||||
if ((iy==0) && (ix>0) ) return 2; // Right atoms
|
||||
if ((iy>0) && (ix==0)) return 1; // Top-middle atoms
|
||||
if ((iy>0) && (ix!=0)) return 3; // Top-right and top-left atoms
|
||||
} else { // iz > 0
|
||||
if((ix==0) && (iy==0)) return 4; // Back atoms
|
||||
if((ix==0) && (iy!=0)) return 5; // Top-back and bottom-back atoms
|
||||
|
|
|
@ -89,7 +89,7 @@ FixManifoldForce::FixManifoldForce(LAMMPS *lmp, int narg, char **arg) :
|
|||
|
||||
double *params = ptr_m->params;
|
||||
for( int i = 0; i < nvars; ++i ){
|
||||
if( was_var( arg[i+4] ) )
|
||||
if (was_var( arg[i+4] ))
|
||||
error->all(FLERR,"Equal-style variables not allowed with fix manifoldforce");
|
||||
|
||||
// Use force->numeric to trigger an error if arg is not a number.
|
||||
|
|
|
@ -84,8 +84,8 @@ FixNVEManifoldRattle::FixNVEManifoldRattle( LAMMPS *lmp, int &narg, char **arg,
|
|||
int error_on_unknown_keyword )
|
||||
: Fix(lmp,narg,arg)
|
||||
{
|
||||
if( lmp->citeme) lmp->citeme->add(cite_fix_nve_manifold_rattle);
|
||||
if( narg < 6 ) error->all(FLERR, "Illegal fix nve/manifold/rattle command");
|
||||
if (lmp->citeme) lmp->citeme->add(cite_fix_nve_manifold_rattle);
|
||||
if (narg < 6 ) error->all(FLERR, "Illegal fix nve/manifold/rattle command");
|
||||
|
||||
// Set all bits/settings:
|
||||
time_integrate = 1;
|
||||
|
@ -135,12 +135,12 @@ FixNVEManifoldRattle::FixNVEManifoldRattle( LAMMPS *lmp, int &narg, char **arg,
|
|||
is_var[i] = 0;
|
||||
}
|
||||
tstrs[i] = new char[len];
|
||||
if( tstrs[i] == NULL ) error->all(FLERR,"Error allocating space for args.");
|
||||
if (tstrs[i] == NULL ) error->all(FLERR,"Error allocating space for args.");
|
||||
strcpy( tstrs[i], arg[i+6] + offset );
|
||||
}
|
||||
|
||||
ptr_m->params = new double[nvars];
|
||||
if( !ptr_m->params ) error->all(FLERR,"Failed to allocate params!");
|
||||
if (!ptr_m->params ) error->all(FLERR,"Failed to allocate params!");
|
||||
for( int i = 0; i < nvars; ++i ){
|
||||
// If param i was variable type, it will be set later...
|
||||
ptr_m->params[i] = is_var[i] ? 0.0 : force->numeric( FLERR, arg[i+6] );
|
||||
|
@ -182,9 +182,9 @@ FixNVEManifoldRattle::~FixNVEManifoldRattle()
|
|||
delete [] tstrs;
|
||||
}
|
||||
|
||||
if( tvars ) delete [] tvars;
|
||||
if( tstyle ) delete [] tstyle;
|
||||
if( is_var ) delete [] is_var;
|
||||
if (tvars ) delete [] tvars;
|
||||
if (tstyle) delete [] tstyle;
|
||||
if (is_var) delete [] is_var;
|
||||
}
|
||||
|
||||
|
||||
|
@ -261,7 +261,7 @@ int FixNVEManifoldRattle::setmask()
|
|||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= FINAL_INTEGRATE;
|
||||
if( nevery > 0 ) mask |= END_OF_STEP;
|
||||
if (nevery > 0) mask |= END_OF_STEP;
|
||||
|
||||
return mask;
|
||||
}
|
||||
|
@ -322,7 +322,7 @@ int FixNVEManifoldRattle::dof(int /*igroup*/)
|
|||
// Make sure that, if there is just no or one atom, no dofs are subtracted,
|
||||
// since for the first atom already 3 dofs are subtracted because of the
|
||||
// centre of mass corrections:
|
||||
if( dofs <= 1 ) dofs = 0;
|
||||
if (dofs <= 1) dofs = 0;
|
||||
stats.dofs_removed = dofs;
|
||||
|
||||
return dofs;
|
||||
|
@ -549,7 +549,7 @@ void FixNVEManifoldRattle::rattle_manifold_x(double *x, double *v,
|
|||
res = infnorm<4>(R);
|
||||
++iters;
|
||||
|
||||
if( (res < tolerance) || (iters >= max_iter) ) break;
|
||||
if ((res < tolerance) || (iters >= max_iter)) break;
|
||||
|
||||
// Update nn and g.
|
||||
gg = ptr_m->g(x);
|
||||
|
|
|
@ -86,7 +86,7 @@ FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
|
|||
{
|
||||
if (lmp->citeme) lmp->citeme->add(cite_fix_nvt_manifold_rattle);
|
||||
|
||||
if( narg < 6 ) error->all(FLERR,"Illegal fix nvt/manifold/rattle command");
|
||||
if (narg < 6 ) error->all(FLERR,"Illegal fix nvt/manifold/rattle command");
|
||||
|
||||
// Set all bits/settings:
|
||||
dof_flag = 1;
|
||||
|
@ -105,7 +105,7 @@ FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
|
|||
while( argi < narg )
|
||||
{
|
||||
if (strcmp( arg[argi], "temp") == 0) {
|
||||
if( argi+3 >= narg )
|
||||
if (argi+3 >= narg)
|
||||
error->all(FLERR,"Keyword 'temp' needs 3 arguments");
|
||||
|
||||
t_start = force->numeric(FLERR, arg[argi+1]);
|
||||
|
@ -116,7 +116,7 @@ FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
|
|||
|
||||
argi += 4;
|
||||
} else if (strcmp( arg[argi], "tchain" ) == 0) {
|
||||
if( argi+1 >= narg )
|
||||
if (argi+1 >= narg)
|
||||
error->all(FLERR,"Keyword 'tchain' needs 1 argument");
|
||||
|
||||
mtchain = force->inumeric(FLERR, arg[argi+1]);
|
||||
|
@ -132,7 +132,7 @@ FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
|
|||
|
||||
reset_dt();
|
||||
|
||||
if( !got_temp ) error->all(FLERR,"Fix nvt/manifold/rattle needs 'temp'!");
|
||||
if (!got_temp ) error->all(FLERR,"Fix nvt/manifold/rattle needs 'temp'!");
|
||||
|
||||
if (t_period < 0.0) {
|
||||
error->all(FLERR,"Fix nvt/manifold/rattle damping parameter must be > 0.0");
|
||||
|
@ -158,7 +158,7 @@ FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
|
|||
"does not exist");
|
||||
}
|
||||
temperature = modify->compute[icompute];
|
||||
if( temperature->tempbias ) which = BIAS;
|
||||
if (temperature->tempbias) which = BIAS;
|
||||
else which = NOBIAS;
|
||||
|
||||
// Set t_freq from t_period
|
||||
|
@ -184,13 +184,13 @@ FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
|
|||
FixNVTManifoldRattle::~FixNVTManifoldRattle()
|
||||
{
|
||||
// Deallocate heap-allocated objects.
|
||||
if( eta ) delete[] eta;
|
||||
if( eta_dot ) delete[] eta_dot;
|
||||
if( eta_dotdot ) delete[] eta_dotdot;
|
||||
if( eta_mass ) delete[] eta_mass;
|
||||
if (eta) delete[] eta;
|
||||
if (eta_dot) delete[] eta_dot;
|
||||
if (eta_dotdot) delete[] eta_dotdot;
|
||||
if (eta_mass) delete[] eta_mass;
|
||||
|
||||
modify->delete_compute(id_temp);
|
||||
if( id_temp ) delete[] id_temp;
|
||||
if (id_temp) delete[] id_temp;
|
||||
}
|
||||
|
||||
|
||||
|
@ -201,7 +201,7 @@ int FixNVTManifoldRattle::setmask()
|
|||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= FINAL_INTEGRATE;
|
||||
if( nevery > 0 ) mask |= END_OF_STEP;
|
||||
if (nevery > 0) mask |= END_OF_STEP;
|
||||
|
||||
return mask;
|
||||
}
|
||||
|
@ -222,7 +222,7 @@ void FixNVTManifoldRattle::init()
|
|||
"does not exist");
|
||||
}
|
||||
temperature = modify->compute[icompute];
|
||||
if( temperature->tempbias ) which = BIAS;
|
||||
if (temperature->tempbias) which = BIAS;
|
||||
else which = NOBIAS;
|
||||
|
||||
}
|
||||
|
@ -341,7 +341,7 @@ void FixNVTManifoldRattle::nh_v_temp()
|
|||
double **v = atom->v;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
if( igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
|
||||
|
||||
|
|
|
@ -42,10 +42,10 @@ public:
|
|||
|
||||
void test()
|
||||
{
|
||||
if( fabs( x(0) - x0 ) > 1e-8 ) err->one(FLERR, "x0 wrong");
|
||||
if( fabs( x(1) - x1 ) > 1e-8 ) err->one(FLERR, "x1 wrong");
|
||||
if( fabs( y(0) - y0 ) > 1e-8 ) err->one(FLERR, "y0 wrong");
|
||||
if( fabs( y(1) - y1 ) > 1e-8 ) err->one(FLERR, "y1 wrong");
|
||||
if (fabs( x(0) - x0 ) > 1e-8 ) err->one(FLERR, "x0 wrong");
|
||||
if (fabs( x(1) - x1 ) > 1e-8 ) err->one(FLERR, "x1 wrong");
|
||||
if (fabs( y(0) - y0 ) > 1e-8 ) err->one(FLERR, "y0 wrong");
|
||||
if (fabs( y(1) - y1 ) > 1e-8 ) err->one(FLERR, "y1 wrong");
|
||||
}
|
||||
|
||||
double get_t_from_x( double xx ) const
|
||||
|
@ -139,8 +139,8 @@ manifold_gaussian_bump::manifold_gaussian_bump(class LAMMPS* lmp,
|
|||
|
||||
manifold_gaussian_bump::~manifold_gaussian_bump()
|
||||
{
|
||||
if( lut_z ) delete lut_z;
|
||||
if( lut_zp ) delete lut_zp;
|
||||
if (lut_z ) delete lut_z;
|
||||
if (lut_zp) delete lut_zp;
|
||||
}
|
||||
|
||||
|
||||
|
@ -349,7 +349,7 @@ void manifold_gaussian_bump::lut_get_z_and_zp( double rr, double &zz,
|
|||
void manifold_gaussian_bump::test_lut()
|
||||
{
|
||||
double x[3], nn[3];
|
||||
if( comm->me != 0 ) return;
|
||||
if (comm->me != 0) return;
|
||||
|
||||
FILE *fp = fopen( "test_lut_gaussian.dat", "w" );
|
||||
double dx = 0.1;
|
||||
|
|
|
@ -123,7 +123,7 @@ thyla_part *manifold_thylakoid::get_thyla_part( const double *x, int * /*err_fla
|
|||
for( std::size_t i = 0; i < parts.size(); ++i ){
|
||||
thyla_part *p = parts[i];
|
||||
if (is_in_domain(p,x)) {
|
||||
if( idx != NULL ) *idx = i;
|
||||
if (idx != NULL) *idx = i;
|
||||
return p;
|
||||
}
|
||||
}
|
||||
|
@ -515,7 +515,7 @@ int manifold_thylakoid::is_in_domain( thyla_part *part, const double *x )
|
|||
(x[1] >= part->ylo) && (x[1] <= part->yhi) &&
|
||||
(x[2] >= part->zlo) && (x[2] <= part->zhi);
|
||||
|
||||
if( !domain_ok ) return false;
|
||||
if (!domain_ok) return false;
|
||||
|
||||
// From here on out, domain is ok.
|
||||
|
||||
|
|
|
@ -38,7 +38,7 @@ thyla_part::thyla_part( int type, double *args, double xlo, double ylo, double z
|
|||
return;
|
||||
}
|
||||
// The others should be 1.
|
||||
if( (args[0] != 1.0) && (args[0] != 0.0) &&
|
||||
if ( (args[0] != 1.0) && (args[0] != 0.0) &&
|
||||
(args[1] != 1.0) && (args[1] != 0.0) &&
|
||||
(args[2] != 1.0) && (args[2] != 0.0) ){
|
||||
err_flag = -1;
|
||||
|
|
|
@ -191,7 +191,7 @@ void PairEDPD::compute(int eflag, int vflag)
|
|||
|
||||
// heat transfer
|
||||
double dQc,dQd,dQr;
|
||||
if( r < cutT[itype][jtype]) {
|
||||
if (r < cutT[itype][jtype]) {
|
||||
double wrT = 1.0 - r/cutT[itype][jtype];
|
||||
wrT = MAX(0.0,MIN(1.0,wrT));
|
||||
wrT = pow(wrT, 0.5*powerT[itype][jtype]);
|
||||
|
|
|
@ -161,7 +161,7 @@ void PairTDPD::compute(int eflag, int vflag)
|
|||
f[i][2] += delz*fpair;
|
||||
|
||||
// chemical concentration transport
|
||||
if( r < cutcc[itype][jtype]) {
|
||||
if (r < cutcc[itype][jtype]) {
|
||||
for(int k=0; k<cc_species; k++) {
|
||||
double wcr = 1.0 - r/cutcc[itype][jtype];
|
||||
wcr = MAX(0,wcr);
|
||||
|
|
|
@ -52,7 +52,7 @@ FixRhok::FixRhok( LAMMPS* inLMP, int inArgc, char** inArgv )
|
|||
if (lmp->citeme) lmp->citeme->add(cite_fix_rhok);
|
||||
|
||||
// Check arguments
|
||||
if( inArgc != 8 )
|
||||
if (inArgc != 8)
|
||||
error->all(FLERR,"Illegal fix rhoKUmbrella command" );
|
||||
|
||||
// Set up fix flags
|
||||
|
@ -104,7 +104,7 @@ void
|
|||
FixRhok::init()
|
||||
{
|
||||
// RESPA boilerplate
|
||||
if( strcmp( update->integrate_style, "respa" ) == 0 )
|
||||
if (strcmp( update->integrate_style, "respa" ) == 0)
|
||||
mNLevelsRESPA = ((Respa *) update->integrate)->nlevels;
|
||||
|
||||
// Count the number of affected particles
|
||||
|
@ -112,7 +112,7 @@ FixRhok::init()
|
|||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
for( int i = 0; i < nlocal; i++ ) { // Iterate through all atoms on this CPU
|
||||
if( mask[i] & groupbit ) { // ...only those affected by this fix
|
||||
if (mask[i] & groupbit) { // ...only those affected by this fix
|
||||
nThisLocal++;
|
||||
}
|
||||
}
|
||||
|
@ -125,7 +125,7 @@ FixRhok::init()
|
|||
void
|
||||
FixRhok::setup( int inVFlag )
|
||||
{
|
||||
if( strcmp( update->integrate_style, "verlet" ) == 0 )
|
||||
if (strcmp( update->integrate_style, "verlet" ) == 0)
|
||||
post_force( inVFlag );
|
||||
else
|
||||
{
|
||||
|
@ -157,7 +157,7 @@ FixRhok::post_force( int /*inVFlag*/ )
|
|||
mRhoKLocal[1] = 0.0;
|
||||
|
||||
for( int i = 0; i < nlocal; i++ ) { // Iterate through all atoms on this CPU
|
||||
if( mask[i] & groupbit ) { // ...only those affected by this fix
|
||||
if (mask[i] & groupbit) { // ...only those affected by this fix
|
||||
|
||||
// rho_k = sum_i exp( - i k.r_i )
|
||||
mRhoKLocal[0] += cos( mK[0]*x[i][0] + mK[1]*x[i][1] + mK[2]*x[i][2] );
|
||||
|
@ -180,7 +180,7 @@ FixRhok::post_force( int /*inVFlag*/ )
|
|||
+ mRhoKGlobal[1]*mRhoKGlobal[1] );
|
||||
|
||||
for( int i = 0; i < nlocal; i++ ) { // Iterate through all atoms on this CPU
|
||||
if( mask[i] & groupbit ) { // ...only those affected by this fix
|
||||
if (mask[i] & groupbit) { // ...only those affected by this fix
|
||||
|
||||
// Calculate forces
|
||||
// U = kappa/2 ( |rho_k| - rho_k^0 )^2
|
||||
|
@ -208,7 +208,7 @@ FixRhok::post_force( int /*inVFlag*/ )
|
|||
void
|
||||
FixRhok::post_force_respa( int inVFlag, int inILevel, int /*inILoop*/ )
|
||||
{
|
||||
if( inILevel == mNLevelsRESPA - 1 )
|
||||
if (inILevel == mNLevelsRESPA - 1)
|
||||
post_force( inVFlag );
|
||||
}
|
||||
|
||||
|
@ -233,7 +233,7 @@ FixRhok::compute_scalar()
|
|||
double
|
||||
FixRhok::compute_vector( int inI )
|
||||
{
|
||||
if( inI == 0 )
|
||||
if (inI == 0)
|
||||
return mRhoKGlobal[0]; // Real part
|
||||
else if( inI == 1 )
|
||||
return mRhoKGlobal[1]; // Imagniary part
|
||||
|
|
|
@ -197,7 +197,7 @@ void PairSRP::compute(int eflag, int vflag)
|
|||
j = jlist[jj];
|
||||
|
||||
// enforce 1-2 exclusions
|
||||
if( (sbmask(j) & exclude) )
|
||||
if ((sbmask(j) & exclude))
|
||||
continue;
|
||||
|
||||
j &= NEIGHMASK;
|
||||
|
@ -274,7 +274,7 @@ void PairSRP::compute(int eflag, int vflag)
|
|||
j = jlist[jj];
|
||||
|
||||
// enforce 1-2 exclusions
|
||||
if( (sbmask(j) & exclude) )
|
||||
if ((sbmask(j) & exclude))
|
||||
continue;
|
||||
|
||||
j &= NEIGHMASK;
|
||||
|
|
|
@ -205,7 +205,7 @@ void PairReaxCOMP::compute(int eflag, int vflag)
|
|||
system->big_box.box_norms[0] = 0;
|
||||
system->big_box.box_norms[1] = 0;
|
||||
system->big_box.box_norms[2] = 0;
|
||||
if( comm->me == 0 ) t_start = MPI_Wtime();
|
||||
if (comm->me == 0 ) t_start = MPI_Wtime();
|
||||
// setup data structures
|
||||
|
||||
setup();
|
||||
|
@ -218,7 +218,7 @@ void PairReaxCOMP::compute(int eflag, int vflag)
|
|||
write_reax_lists();
|
||||
|
||||
// timing for filling in the reax lists
|
||||
if( comm->me == 0 ) {
|
||||
if (comm->me == 0) {
|
||||
t_end = MPI_Wtime();
|
||||
data->timing.nbrs = t_end - t_start;
|
||||
}
|
||||
|
|
|
@ -113,7 +113,7 @@ void Add_dBond_to_ForcesOMP( reax_system *system, int i, int pj,
|
|||
|
||||
rvec_Add(workspace->forceReduction[reductionOffset+i],temp );
|
||||
|
||||
if( system->pair_ptr->vflag_atom) {
|
||||
if (system->pair_ptr->vflag_atom) {
|
||||
rvec_Scale(fi_tmp, -1.0, temp);
|
||||
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,-1., system->my_atoms[j].x );
|
||||
|
||||
|
@ -147,7 +147,7 @@ void Add_dBond_to_ForcesOMP( reax_system *system, int i, int pj,
|
|||
|
||||
rvec_Add(workspace->forceReduction[reductionOffset+j],temp );
|
||||
|
||||
if( system->pair_ptr->vflag_atom) {
|
||||
if (system->pair_ptr->vflag_atom) {
|
||||
rvec_Scale(fj_tmp, -1.0, temp);
|
||||
rvec_ScaledSum( delji, 1., system->my_atoms[j].x,-1., system->my_atoms[i].x );
|
||||
|
||||
|
@ -171,7 +171,7 @@ void Add_dBond_to_ForcesOMP( reax_system *system, int i, int pj,
|
|||
|
||||
rvec_Add(workspace->forceReduction[reductionOffset+k],temp );
|
||||
|
||||
if( system->pair_ptr->vflag_atom ) {
|
||||
if (system->pair_ptr->vflag_atom) {
|
||||
rvec_Scale(fk_tmp, -1.0, temp);
|
||||
rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x);
|
||||
|
||||
|
@ -201,7 +201,7 @@ void Add_dBond_to_ForcesOMP( reax_system *system, int i, int pj,
|
|||
|
||||
rvec_Add(workspace->forceReduction[reductionOffset+k],temp );
|
||||
|
||||
if( system->pair_ptr->vflag_atom ) {
|
||||
if (system->pair_ptr->vflag_atom) {
|
||||
rvec_Scale(fk_tmp, -1.0, temp);
|
||||
rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x);
|
||||
|
||||
|
@ -313,7 +313,7 @@ void Add_dBond_to_Forces_NPTOMP( reax_system *system, int i, int pj,
|
|||
rvec_Add(workspace->forceReduction[reductionOffset+k],temp );
|
||||
|
||||
/* pressure */
|
||||
if( k != i ) {
|
||||
if (k != i) {
|
||||
ivec_Sum( rel_box, nbr_k->rel_box, nbr_j->rel_box ); //rel_box(k, i)
|
||||
rvec_iMultiply( ext_press, rel_box, temp );
|
||||
rvec_Add( workspace->my_ext_pressReduction[tid], ext_press );
|
||||
|
@ -490,10 +490,10 @@ void BOOMP( reax_system *system, control_params * /* control */, simulation_data
|
|||
if(type_j < 0) continue;
|
||||
bo_ij = &( bonds->select.bond_list[pj].bo_data );
|
||||
|
||||
if( i < j || workspace->bond_mark[j] > 3) {
|
||||
if (i < j || workspace->bond_mark[j] > 3) {
|
||||
twbp = &( system->reax_param.tbp[type_i][type_j] );
|
||||
|
||||
if( twbp->ovc < 0.001 && twbp->v13cor < 0.001 ) {
|
||||
if (twbp->ovc < 0.001 && twbp->v13cor < 0.001) {
|
||||
bo_ij->C1dbo = 1.000000;
|
||||
bo_ij->C2dbo = 0.000000;
|
||||
bo_ij->C3dbo = 0.000000;
|
||||
|
@ -515,7 +515,7 @@ void BOOMP( reax_system *system, control_params * /* control */, simulation_data
|
|||
Deltap_boc_j = workspace->Deltap_boc[j];
|
||||
|
||||
/* on page 1 */
|
||||
if( twbp->ovc >= 0.001 ) {
|
||||
if (twbp->ovc >= 0.001) {
|
||||
/* Correction for overcoordination */
|
||||
exp_p1i = exp( -p_boc1 * Deltap_i );
|
||||
exp_p2i = exp( -p_boc2 * Deltap_i );
|
||||
|
@ -556,7 +556,7 @@ void BOOMP( reax_system *system, control_params * /* control */, simulation_data
|
|||
Cf1_ij = Cf1_ji = 0.0;
|
||||
}
|
||||
|
||||
if( twbp->v13cor >= 0.001 ) {
|
||||
if (twbp->v13cor >= 0.001) {
|
||||
/* Correction for 1-3 bond orders */
|
||||
exp_f4 =exp(-(twbp->p_boc4 * SQR( bo_ij->BO ) -
|
||||
Deltap_boc_i) * twbp->p_boc3 + twbp->p_boc5);
|
||||
|
@ -607,13 +607,13 @@ void BOOMP( reax_system *system, control_params * /* control */, simulation_data
|
|||
}
|
||||
|
||||
/* neglect bonds that are < 1e-10 */
|
||||
if( bo_ij->BO < 1e-10 )
|
||||
if (bo_ij->BO < 1e-10)
|
||||
bo_ij->BO = 0.0;
|
||||
if( bo_ij->BO_s < 1e-10 )
|
||||
if (bo_ij->BO_s < 1e-10)
|
||||
bo_ij->BO_s = 0.0;
|
||||
if( bo_ij->BO_pi < 1e-10 )
|
||||
if (bo_ij->BO_pi < 1e-10)
|
||||
bo_ij->BO_pi = 0.0;
|
||||
if( bo_ij->BO_pi2 < 1e-10 )
|
||||
if (bo_ij->BO_pi2 < 1e-10)
|
||||
bo_ij->BO_pi2 = 0.0;
|
||||
|
||||
workspace->total_bond_order[i] += bo_ij->BO; //now keeps total_BO
|
||||
|
@ -654,7 +654,7 @@ void BOOMP( reax_system *system, control_params * /* control */, simulation_data
|
|||
type_j = system->my_atoms[j].type;
|
||||
if(type_j < 0) continue;
|
||||
|
||||
if( i < j || workspace->bond_mark[j] > 3) {
|
||||
if (i < j || workspace->bond_mark[j] > 3) {
|
||||
// Computed in previous for-loop
|
||||
} else {
|
||||
/* We only need to update bond orders from bo_ji
|
||||
|
@ -705,7 +705,7 @@ void BOOMP( reax_system *system, control_params * /* control */, simulation_data
|
|||
workspace->Clp[j] = 2.0 * p_lp1 * explp1 * (2.0 + workspace->vlpex[j]);
|
||||
workspace->dDelta_lp[j] = workspace->Clp[j];
|
||||
|
||||
if( sbp_j->mass > 21.0 ) {
|
||||
if (sbp_j->mass > 21.0) {
|
||||
workspace->nlp_temp[j] = 0.5 * (sbp_j->valency_e - sbp_j->valency);
|
||||
workspace->Delta_lp_temp[j] = sbp_j->nlp_opt - workspace->nlp_temp[j];
|
||||
workspace->dDelta_lp_temp[j] = 0.;
|
||||
|
|
|
@ -99,9 +99,9 @@ void BondsOMP( reax_system *system, control_params * /* control */,
|
|||
for (pj = start_i; pj < end_i; ++pj) {
|
||||
j = bonds->select.bond_list[pj].nbr;
|
||||
|
||||
if( system->my_atoms[i].orig_id > system->my_atoms[j].orig_id ) continue;
|
||||
if (system->my_atoms[i].orig_id > system->my_atoms[j].orig_id) continue;
|
||||
|
||||
if( system->my_atoms[i].orig_id == system->my_atoms[j].orig_id ) {
|
||||
if (system->my_atoms[i].orig_id == system->my_atoms[j].orig_id) {
|
||||
if (system->my_atoms[j].x[2] < system->my_atoms[i].x[2]) continue;
|
||||
if (system->my_atoms[j].x[2] == system->my_atoms[i].x[2] &&
|
||||
system->my_atoms[j].x[1] < system->my_atoms[i].x[1]) continue;
|
||||
|
|
|
@ -64,7 +64,7 @@ void Init_Force_FunctionsOMP( control_params *control )
|
|||
Interaction_Functions[2] = Atom_EnergyOMP; //Dummy_Interaction;
|
||||
Interaction_Functions[3] = Valence_AnglesOMP; //Dummy_Interaction;
|
||||
Interaction_Functions[4] = Torsion_AnglesOMP; //Dummy_Interaction;
|
||||
if( control->hbond_cut > 0 )
|
||||
if (control->hbond_cut > 0)
|
||||
Interaction_Functions[5] = Hydrogen_BondsOMP;
|
||||
else Interaction_Functions[5] = Dummy_Interaction;
|
||||
Interaction_Functions[6] = Dummy_Interaction; //empty
|
||||
|
@ -113,7 +113,7 @@ void Compute_NonBonded_ForcesOMP( reax_system *system, control_params *control,
|
|||
startTimeBase = MPI_Wtime();
|
||||
#endif
|
||||
|
||||
if( control->tabulate == 0 )
|
||||
if (control->tabulate == 0)
|
||||
vdW_Coulomb_Energy_OMP( system, control, data, workspace,
|
||||
lists, out_control );
|
||||
else
|
||||
|
@ -275,7 +275,7 @@ void Validate_ListsOMP( reax_system *system, storage * /*workspace */, reax_list
|
|||
{
|
||||
|
||||
/* bond list */
|
||||
if( N > 0 ) {
|
||||
if (N > 0) {
|
||||
bonds = *lists + BONDS;
|
||||
|
||||
#if defined(_OPENMP)
|
||||
|
@ -284,11 +284,11 @@ void Validate_ListsOMP( reax_system *system, storage * /*workspace */, reax_list
|
|||
for( i = 0; i < N; ++i ) {
|
||||
system->my_atoms[i].num_bonds = MAX(Num_Entries(i,bonds)*2, MIN_BONDS);
|
||||
|
||||
if( i < N-1 )
|
||||
if (i < N-1)
|
||||
comp = Start_Index(i+1, bonds);
|
||||
else comp = bonds->num_intrs;
|
||||
|
||||
if( End_Index(i, bonds) > comp ) {
|
||||
if (End_Index(i, bonds) > comp) {
|
||||
fprintf( stderr, "step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
|
||||
step, i, End_Index(i,bonds), comp );
|
||||
MPI_Abort( comm, INSUFFICIENT_MEMORY );
|
||||
|
@ -298,7 +298,7 @@ void Validate_ListsOMP( reax_system *system, storage * /*workspace */, reax_list
|
|||
|
||||
|
||||
/* hbonds list */
|
||||
if( numH > 0 ) {
|
||||
if (numH > 0) {
|
||||
hbonds = *lists + HBONDS;
|
||||
|
||||
#if defined(_OPENMP)
|
||||
|
@ -306,15 +306,15 @@ void Validate_ListsOMP( reax_system *system, storage * /*workspace */, reax_list
|
|||
#endif
|
||||
for( i = 0; i < n; ++i ) {
|
||||
Hindex = system->my_atoms[i].Hindex;
|
||||
if( Hindex > -1 ) {
|
||||
if (Hindex > -1) {
|
||||
system->my_atoms[i].num_hbonds =
|
||||
(int)(MAX( Num_Entries(Hindex, hbonds)*saferzone, MIN_HBONDS ));
|
||||
|
||||
if( Hindex < numH-1 )
|
||||
if (Hindex < numH-1)
|
||||
comp = Start_Index(Hindex+1, hbonds);
|
||||
else comp = hbonds->num_intrs;
|
||||
|
||||
if( End_Index(Hindex, hbonds) > comp ) {
|
||||
if (End_Index(Hindex, hbonds) > comp) {
|
||||
fprintf(stderr,"step%d-hbondchk failed: H=%d end(H)=%d str(H+1)=%d\n",
|
||||
step, Hindex, End_Index(Hindex,hbonds), comp );
|
||||
MPI_Abort( comm, INSUFFICIENT_MEMORY );
|
||||
|
@ -402,7 +402,7 @@ void Init_Forces_noQEq_OMP( reax_system *system, control_params *control,
|
|||
// #pragma omp critical
|
||||
// {
|
||||
// btop_i = End_Index(i, bonds);
|
||||
// if( BOp(workspace, bonds, control->bo_cut, i, btop_i, nbr_pj, sbp_i, sbp_j, twbp) ) {
|
||||
// if (BOp(workspace, bonds, control->bo_cut, i, btop_i, nbr_pj, sbp_i, sbp_j, twbp)) {
|
||||
// num_bonds++;
|
||||
// btop_i++;
|
||||
// Set_End_Index(i, btop_i, bonds);
|
||||
|
@ -418,19 +418,19 @@ void Init_Forces_noQEq_OMP( reax_system *system, control_params *control,
|
|||
double BO, BO_s, BO_pi, BO_pi2;
|
||||
double bo_cut = control->bo_cut;
|
||||
|
||||
if( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0 ) {
|
||||
if (sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0) {
|
||||
C12 = twbp->p_bo1 * pow( nbr_pj->d / twbp->r_s, twbp->p_bo2 );
|
||||
BO_s = (1.0 + bo_cut) * exp( C12 );
|
||||
}
|
||||
else BO_s = C12 = 0.0;
|
||||
|
||||
if( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 ) {
|
||||
if (sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0) {
|
||||
C34 = twbp->p_bo3 * pow( nbr_pj->d / twbp->r_p, twbp->p_bo4 );
|
||||
BO_pi = exp( C34 );
|
||||
}
|
||||
else BO_pi = C34 = 0.0;
|
||||
|
||||
if( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 ) {
|
||||
if (sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0) {
|
||||
C56 = twbp->p_bo5 * pow( nbr_pj->d / twbp->r_pp, twbp->p_bo6 );
|
||||
BO_pi2= exp( C56 );
|
||||
}
|
||||
|
|
|
@ -110,7 +110,7 @@ void Hydrogen_BondsOMP( reax_system *system, control_params *control,
|
|||
// for( j = 0; j < system->n; ++j )
|
||||
for( j = ifrom; j < ito; ++j ) {
|
||||
/* j has to be of type H */
|
||||
if( system->reax_param.sbp[system->my_atoms[j].type].p_hbond == 1 ) {
|
||||
if (system->reax_param.sbp[system->my_atoms[j].type].p_hbond == 1) {
|
||||
/*set j's variables */
|
||||
type_j = system->my_atoms[j].type;
|
||||
start_j = Start_Index(j, bonds);
|
||||
|
@ -146,7 +146,7 @@ void Hydrogen_BondsOMP( reax_system *system, control_params *control,
|
|||
pbond_ij = &( bonds->select.bond_list[pi] );
|
||||
i = pbond_ij->nbr;
|
||||
|
||||
if( system->my_atoms[i].orig_id != system->my_atoms[k].orig_id ) {
|
||||
if (system->my_atoms[i].orig_id != system->my_atoms[k].orig_id) {
|
||||
bo_ij = &(pbond_ij->bo_data);
|
||||
type_i = system->my_atoms[i].type;
|
||||
if(type_i < 0) continue;
|
||||
|
@ -179,7 +179,7 @@ void Hydrogen_BondsOMP( reax_system *system, control_params *control,
|
|||
/* hydrogen bond forces */
|
||||
bo_ij->Cdbo += CEhb1; // dbo term
|
||||
|
||||
if( control->virial == 0 ) {
|
||||
if (control->virial == 0) {
|
||||
// dcos terms
|
||||
rvec_ScaledAdd(workspace->forceReduction[reductionOffset+i], +CEhb2, dcos_theta_di );
|
||||
rvec_ScaledAdd(workspace->forceReduction[reductionOffset+j], +CEhb2, dcos_theta_dj );
|
||||
|
|
|
@ -65,7 +65,7 @@ int Init_ListsOMP( reax_system *system, control_params *control,
|
|||
Estimate_Storages( system, control, lists,
|
||||
&Htop, hb_top, bond_top, &num_3body, comm );
|
||||
|
||||
if( control->hbond_cut > 0 ) {
|
||||
if (control->hbond_cut > 0) {
|
||||
/* init H indexes */
|
||||
total_hbonds = 0;
|
||||
for( i = 0; i < system->n; ++i ) {
|
||||
|
@ -127,7 +127,7 @@ void InitializeOMP( reax_system *system, control_params *control,
|
|||
char msg[MAX_STR];
|
||||
|
||||
|
||||
if( Init_MPI_Datatypes(system, workspace, mpi_data, comm, msg) == FAILURE ) {
|
||||
if (Init_MPI_Datatypes(system, workspace, mpi_data, comm, msg) == FAILURE) {
|
||||
fprintf( stderr, "p%d: init_mpi_datatypes: could not create datatypes\n",
|
||||
system->my_rank );
|
||||
fprintf( stderr, "p%d: mpi_data couldn't be initialized! terminating.\n",
|
||||
|
@ -142,15 +142,15 @@ void InitializeOMP( reax_system *system, control_params *control,
|
|||
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
|
||||
}
|
||||
|
||||
if( Init_Simulation_Data( system, control, data, msg ) == FAILURE ) {
|
||||
if (Init_Simulation_Data( system, control, data, msg ) == FAILURE) {
|
||||
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
|
||||
fprintf( stderr, "p%d: sim_data couldn't be initialized! terminating.\n",
|
||||
system->my_rank );
|
||||
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
|
||||
}
|
||||
|
||||
if( Init_Workspace( system, control, workspace, mpi_data->world, msg ) ==
|
||||
FAILURE ) {
|
||||
if (Init_Workspace( system, control, workspace, mpi_data->world, msg ) ==
|
||||
FAILURE) {
|
||||
fprintf( stderr, "p%d:init_workspace: not enough memory\n",
|
||||
system->my_rank );
|
||||
fprintf( stderr, "p%d:workspace couldn't be initialized! terminating.\n",
|
||||
|
@ -158,23 +158,23 @@ void InitializeOMP( reax_system *system, control_params *control,
|
|||
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
|
||||
}
|
||||
|
||||
if( Init_ListsOMP( system, control, data, workspace, lists, mpi_data, msg ) ==
|
||||
FAILURE ) {
|
||||
if (Init_ListsOMP( system, control, data, workspace, lists, mpi_data, msg ) ==
|
||||
FAILURE) {
|
||||
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
|
||||
fprintf( stderr, "p%d: system could not be initialized! terminating.\n",
|
||||
system->my_rank );
|
||||
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
|
||||
}
|
||||
|
||||
if( Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE) {
|
||||
if (Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE) {
|
||||
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
|
||||
fprintf( stderr, "p%d: could not open output files! terminating...\n",
|
||||
system->my_rank );
|
||||
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
|
||||
}
|
||||
|
||||
if( control->tabulate ) {
|
||||
if( Init_Lookup_Tables( system, control, workspace, mpi_data, msg ) == FAILURE ) {
|
||||
if (control->tabulate) {
|
||||
if (Init_Lookup_Tables( system, control, workspace, mpi_data, msg ) == FAILURE) {
|
||||
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
|
||||
fprintf( stderr, "p%d: couldn't create lookup table! terminating.\n",
|
||||
system->my_rank );
|
||||
|
|
|
@ -128,24 +128,24 @@ void Atom_EnergyOMP( reax_system *system, control_params * /* control */,
|
|||
if(numbonds > 0) workspace->CdDelta[i] += CElp; // lp - 1st term
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag)
|
||||
if (system->pair_ptr->evflag)
|
||||
pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, i, system->n, 1,
|
||||
e_lp, 0.0, 0.0, 0.0, 0.0, 0.0, thr);
|
||||
|
||||
/* correction for C2 */
|
||||
if( p_lp3 > 0.001 && !strcmp(system->reax_param.sbp[type_i].name, "C") )
|
||||
if (p_lp3 > 0.001 && !strcmp(system->reax_param.sbp[type_i].name, "C"))
|
||||
for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj ) {
|
||||
j = bonds->select.bond_list[pj].nbr;
|
||||
type_j = system->my_atoms[j].type;
|
||||
if(type_j < 0) continue;
|
||||
|
||||
if( !strcmp( system->reax_param.sbp[type_j].name, "C" ) ) {
|
||||
if (!strcmp( system->reax_param.sbp[type_j].name, "C" )) {
|
||||
twbp = &( system->reax_param.tbp[type_i][type_j]);
|
||||
bo_ij = &( bonds->select.bond_list[pj].bo_data );
|
||||
Di = workspace->Delta[i];
|
||||
vov3 = bo_ij->BO - Di - 0.040*pow(Di, 4.);
|
||||
|
||||
if( vov3 > 3. ) {
|
||||
if (vov3 > 3.) {
|
||||
total_Elp += e_lph = p_lp3 * SQR(vov3-3.0);
|
||||
|
||||
deahu2dbo = 2.*p_lp3*(vov3 - 3.);
|
||||
|
@ -155,7 +155,7 @@ void Atom_EnergyOMP( reax_system *system, control_params * /* control */,
|
|||
workspace->CdDelta[i] += deahu2dsbo;
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag)
|
||||
if (system->pair_ptr->evflag)
|
||||
pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, i, j, system->n, 1,
|
||||
e_lph, 0.0, 0.0, 0.0, 0.0, 0.0, thr);
|
||||
}
|
||||
|
@ -172,7 +172,7 @@ void Atom_EnergyOMP( reax_system *system, control_params * /* control */,
|
|||
sbp_i = &(system->reax_param.sbp[ type_i ]);
|
||||
|
||||
/* over-coordination energy */
|
||||
if( sbp_i->mass > 21.0 )
|
||||
if (sbp_i->mass > 21.0)
|
||||
dfvl = 0.0;
|
||||
else dfvl = 1.0; // only for 1st-row elements
|
||||
|
||||
|
|
|
@ -208,7 +208,7 @@ void vdW_Coulomb_Energy_OMP( reax_system *system, control_params *control,
|
|||
( dTap - Tap * r_ij / dr3gamij_1 ) / dr3gamij_3;
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
|
||||
if (system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
|
||||
pe_vdw = Tap * (e_vdW + e_core + e_lg);
|
||||
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
|
||||
-1., system->my_atoms[j].x );
|
||||
|
@ -218,7 +218,7 @@ void vdW_Coulomb_Energy_OMP( reax_system *system, control_params *control,
|
|||
delij[0], delij[1], delij[2], thr);
|
||||
}
|
||||
|
||||
if( control->virial == 0 ) {
|
||||
if (control->virial == 0) {
|
||||
rvec_ScaledAdd( workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec );
|
||||
rvec_ScaledAdd( workspace->forceReduction[reductionOffset+j],
|
||||
+(CEvd + CEclmb), nbr_pj->dvec );
|
||||
|
@ -327,7 +327,7 @@ void Tabulated_vdW_Coulomb_Energy_OMP(reax_system *system,control_params *contro
|
|||
|
||||
/* Cubic Spline Interpolation */
|
||||
r = (int)(r_ij * t->inv_dx);
|
||||
if( r == 0 ) ++r;
|
||||
if (r == 0) ++r;
|
||||
base = (double)(r+1) * t->dx;
|
||||
dif = r_ij - base;
|
||||
|
||||
|
@ -349,7 +349,7 @@ void Tabulated_vdW_Coulomb_Energy_OMP(reax_system *system,control_params *contro
|
|||
CEclmb *= system->my_atoms[i].q * system->my_atoms[j].q;
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
|
||||
if (system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
|
||||
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
|
||||
-1., system->my_atoms[j].x );
|
||||
f_tmp = -(CEvd + CEclmb);
|
||||
|
@ -357,7 +357,7 @@ void Tabulated_vdW_Coulomb_Energy_OMP(reax_system *system,control_params *contro
|
|||
f_tmp, delij[0], delij[1], delij[2], thr);
|
||||
}
|
||||
|
||||
if( control->virial == 0 ) {
|
||||
if (control->virial == 0) {
|
||||
rvec_ScaledAdd( workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec );
|
||||
rvec_ScaledAdd( workspace->forceReduction[froffset+j],
|
||||
+(CEvd + CEclmb), nbr_pj->dvec );
|
||||
|
|
|
@ -198,7 +198,7 @@ void Torsion_AnglesOMP( reax_system *system, control_params *control,
|
|||
sin_ijk = sin( theta_ijk );
|
||||
cos_ijk = cos( theta_ijk );
|
||||
//tan_ijk_i = 1. / tan( theta_ijk );
|
||||
if( sin_ijk >= 0 && sin_ijk <= MIN_SINE )
|
||||
if (sin_ijk >= 0 && sin_ijk <= MIN_SINE)
|
||||
tan_ijk_i = cos_ijk / MIN_SINE;
|
||||
else if( sin_ijk <= 0 && sin_ijk >= -MIN_SINE )
|
||||
tan_ijk_i = cos_ijk / -MIN_SINE;
|
||||
|
@ -237,7 +237,7 @@ void Torsion_AnglesOMP( reax_system *system, control_params *control,
|
|||
sin_jkl = sin( theta_jkl );
|
||||
cos_jkl = cos( theta_jkl );
|
||||
//tan_jkl_i = 1. / tan( theta_jkl );
|
||||
if( sin_jkl >= 0 && sin_jkl <= MIN_SINE )
|
||||
if (sin_jkl >= 0 && sin_jkl <= MIN_SINE)
|
||||
tan_jkl_i = cos_jkl / MIN_SINE;
|
||||
else if( sin_jkl <= 0 && sin_jkl >= -MIN_SINE )
|
||||
tan_jkl_i = cos_jkl / -MIN_SINE;
|
||||
|
@ -338,7 +338,7 @@ void Torsion_AnglesOMP( reax_system *system, control_params *control,
|
|||
//bo_kl->Cdbo += (CEtors6 + CEconj3);
|
||||
bo_kl->CdboReduction[tid] += (CEtors6 + CEconj3);
|
||||
|
||||
if( control->virial == 0 ) {
|
||||
if (control->virial == 0) {
|
||||
/* dcos_theta_ijk */
|
||||
rvec_ScaledAdd( workspace->f[j],
|
||||
CEtors7 + CEconj4, p_ijk->dcos_dj );
|
||||
|
@ -417,7 +417,7 @@ void Torsion_AnglesOMP( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* tally into per-atom virials */
|
||||
if( system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
|
||||
if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
|
||||
|
||||
// acquire vectors
|
||||
rvec_ScaledSum( delil, 1., system->my_atoms[l].x,
|
||||
|
|
|
@ -319,13 +319,13 @@ void Valence_AnglesOMP( reax_system *system, control_params *control,
|
|||
SBO = SBOp + (1 - prod_SBO) * (-workspace->Delta_boc[j] - p_val8 * vlpadj);
|
||||
dSBO1 = -8 * prod_SBO * ( workspace->Delta_boc[j] + p_val8 * vlpadj );
|
||||
|
||||
if( SBO <= 0 )
|
||||
if (SBO <= 0)
|
||||
SBO2 = 0, CSBO2 = 0;
|
||||
else if( SBO > 0 && SBO <= 1 ) {
|
||||
else if (SBO > 0 && SBO <= 1) {
|
||||
SBO2 = pow( SBO, p_val9 );
|
||||
CSBO2 = p_val9 * pow( SBO, p_val9 - 1 );
|
||||
}
|
||||
else if( SBO > 1 && SBO < 2 ) {
|
||||
else if (SBO > 1 && SBO < 2) {
|
||||
SBO2 = 2 - pow( 2-SBO, p_val9 );
|
||||
CSBO2 = p_val9 * pow( 2 - SBO, p_val9 - 1 );
|
||||
}
|
||||
|
@ -398,7 +398,7 @@ void Valence_AnglesOMP( reax_system *system, control_params *control,
|
|||
p_ijk->theta = theta;
|
||||
|
||||
sin_theta = sin( theta );
|
||||
if( sin_theta < 1.0e-5 )
|
||||
if (sin_theta < 1.0e-5)
|
||||
sin_theta = 1.0e-5;
|
||||
|
||||
++my_offset; // add this to the list of 3-body interactions
|
||||
|
@ -412,7 +412,7 @@ void Valence_AnglesOMP( reax_system *system, control_params *control,
|
|||
|
||||
for (cnt = 0; cnt < thbh->cnt; ++cnt) {
|
||||
|
||||
if( fabs(thbh->prm[cnt].p_val1) > 0.001 ) {
|
||||
if (fabs(thbh->prm[cnt].p_val1) > 0.001) {
|
||||
thbp = &( thbh->prm[cnt] );
|
||||
|
||||
/* ANGLE ENERGY */
|
||||
|
@ -536,7 +536,7 @@ void Valence_AnglesOMP( reax_system *system, control_params *control,
|
|||
bo_jt->Cdbopi2 += CEval5;
|
||||
}
|
||||
|
||||
if( control->virial == 0 ) {
|
||||
if (control->virial == 0) {
|
||||
rvec_ScaledAdd( workspace->f[j], CEval8, p_ijk->dcos_dj );
|
||||
rvec_ScaledAdd( workspace->forceReduction[reductionOffset+i],
|
||||
CEval8, p_ijk->dcos_di );
|
||||
|
@ -561,7 +561,7 @@ void Valence_AnglesOMP( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* tally into per-atom virials */
|
||||
if( system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
|
||||
if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
|
||||
|
||||
/* Acquire vectors */
|
||||
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
|
||||
|
@ -575,10 +575,10 @@ void Valence_AnglesOMP( reax_system *system, control_params *control,
|
|||
|
||||
eng_tmp = e_ang + e_pen + e_coa;
|
||||
|
||||
if( system->pair_ptr->evflag)
|
||||
if (system->pair_ptr->evflag)
|
||||
pair_reax_ptr->ev_tally_thr_proxy(system->pair_ptr, j, j, system->N, 1,
|
||||
eng_tmp, 0.0, 0.0, 0.0, 0.0, 0.0, thr);
|
||||
if( system->pair_ptr->vflag_atom)
|
||||
if (system->pair_ptr->vflag_atom)
|
||||
// NEED TO MAKE AN OMP VERSION OF THIS CALL!
|
||||
system->pair_ptr->v_tally3( i, j, k, fi_tmp, fk_tmp, delij, delkj);
|
||||
}
|
||||
|
@ -598,9 +598,9 @@ void Valence_AnglesOMP( reax_system *system, control_params *control,
|
|||
data->my_en.e_pen = total_Epen;
|
||||
data->my_en.e_coa = total_Ecoa;
|
||||
|
||||
if( num_thb_intrs >= thb_intrs->num_intrs * DANGER_ZONE ) {
|
||||
if (num_thb_intrs >= thb_intrs->num_intrs * DANGER_ZONE) {
|
||||
workspace->realloc.num_3body = num_thb_intrs * TWICE;
|
||||
if( num_thb_intrs > thb_intrs->num_intrs ) {
|
||||
if (num_thb_intrs > thb_intrs->num_intrs) {
|
||||
fprintf( stderr, "step%d-ran out of space on angle_list: top=%d, max=%d",
|
||||
data->step, num_thb_intrs, thb_intrs->num_intrs );
|
||||
MPI_Abort( MPI_COMM_WORLD, INSUFFICIENT_MEMORY );
|
||||
|
|
|
@ -252,7 +252,7 @@ void PairQUIP::coeff(int narg, char **arg)
|
|||
}
|
||||
}
|
||||
|
||||
if( narg != (4+n) ) {
|
||||
if (narg != (4+n)) {
|
||||
char str[1024];
|
||||
sprintf(str,"Number of arguments %d is not correct, it should be %d", narg, 4+n);
|
||||
error->all(FLERR,str);
|
||||
|
@ -271,7 +271,7 @@ void PairQUIP::coeff(int narg, char **arg)
|
|||
|
||||
for (int i = 4; i < narg; i++) {
|
||||
|
||||
if( 0 == sscanf(arg[i],"%d",&map[i-4])) {
|
||||
if (0 == sscanf(arg[i],"%d",&map[i-4])) {
|
||||
char str[1024];
|
||||
sprintf(str,"Incorrect atomic number %s at position %d",arg[i],i);
|
||||
error->all(FLERR,str);
|
||||
|
|
|
@ -137,9 +137,9 @@ PairReaxC::~PairReaxC()
|
|||
|
||||
// deallocate reax data-structures
|
||||
|
||||
if( control->tabulate ) Deallocate_Lookup_Tables( system );
|
||||
if (control->tabulate ) Deallocate_Lookup_Tables( system);
|
||||
|
||||
if( control->hbond_cut > 0 ) Delete_List( lists+HBONDS, world );
|
||||
if (control->hbond_cut > 0 ) Delete_List( lists+HBONDS, world);
|
||||
Delete_List( lists+BONDS, world );
|
||||
Delete_List( lists+THREE_BODIES, world );
|
||||
Delete_List( lists+FAR_NBRS, world );
|
||||
|
@ -157,7 +157,7 @@ PairReaxC::~PairReaxC()
|
|||
memory->destroy( mpi_data );
|
||||
|
||||
// deallocate interface storage
|
||||
if( allocated ) {
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
memory->destroy(cutghost);
|
||||
|
@ -510,7 +510,7 @@ void PairReaxC::compute(int eflag, int vflag)
|
|||
system->big_box.box_norms[0] = 0;
|
||||
system->big_box.box_norms[1] = 0;
|
||||
system->big_box.box_norms[2] = 0;
|
||||
if( comm->me == 0 ) t_start = MPI_Wtime();
|
||||
if (comm->me == 0 ) t_start = MPI_Wtime();
|
||||
|
||||
// setup data structures
|
||||
|
||||
|
@ -519,7 +519,7 @@ void PairReaxC::compute(int eflag, int vflag)
|
|||
Reset( system, control, data, workspace, &lists, world );
|
||||
workspace->realloc.num_far = write_reax_lists();
|
||||
// timing for filling in the reax lists
|
||||
if( comm->me == 0 ) {
|
||||
if (comm->me == 0) {
|
||||
t_end = MPI_Wtime();
|
||||
data->timing.nbrs = t_end - t_start;
|
||||
}
|
||||
|
@ -683,7 +683,7 @@ int PairReaxC::estimate_reax_lists()
|
|||
j &= NEIGHMASK;
|
||||
get_distance( x[j], x[i], &d_sqr, &dvec );
|
||||
|
||||
if( d_sqr <= SQR(control->nonb_cut) )
|
||||
if (d_sqr <= SQR(control->nonb_cut))
|
||||
++num_nbrs;
|
||||
}
|
||||
}
|
||||
|
|
|
@ -122,7 +122,7 @@ void DeAllocate_Workspace( control_params * /*control*/, storage *workspace )
|
|||
{
|
||||
int i;
|
||||
|
||||
if( !workspace->allocated )
|
||||
if (!workspace->allocated)
|
||||
return;
|
||||
|
||||
workspace->allocated = 0;
|
||||
|
@ -338,13 +338,13 @@ static int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds,
|
|||
|
||||
total_hbonds = 0;
|
||||
for( i = 0; i < system->n; ++i )
|
||||
if( (system->my_atoms[i].Hindex) >= 0 ) {
|
||||
if ((system->my_atoms[i].Hindex) >= 0) {
|
||||
total_hbonds += system->my_atoms[i].num_hbonds;
|
||||
}
|
||||
total_hbonds = (int)(MAX( total_hbonds*saferzone, mincap*MIN_HBONDS ));
|
||||
|
||||
Delete_List( hbonds, comm );
|
||||
if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND, hbonds, comm ) ) {
|
||||
if (!Make_List( system->Hcap, total_hbonds, TYP_HBOND, hbonds, comm )) {
|
||||
fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
|
||||
MPI_Abort( comm, INSUFFICIENT_MEMORY );
|
||||
}
|
||||
|
@ -429,10 +429,10 @@ void ReAllocate( reax_system *system, control_params *control,
|
|||
system->total_cap = MAX( (int)(system->N * safezone), mincap );
|
||||
}
|
||||
|
||||
if( Nflag ) {
|
||||
if (Nflag) {
|
||||
/* system */
|
||||
ret = Allocate_System( system, system->local_cap, system->total_cap, msg );
|
||||
if( ret != SUCCESS ) {
|
||||
if (ret != SUCCESS) {
|
||||
fprintf( stderr, "not enough space for atom_list: total_cap=%d",
|
||||
system->total_cap );
|
||||
fprintf( stderr, "terminating...\n" );
|
||||
|
@ -443,7 +443,7 @@ void ReAllocate( reax_system *system, control_params *control,
|
|||
DeAllocate_Workspace( control, workspace );
|
||||
ret = Allocate_Workspace( system, control, workspace, system->local_cap,
|
||||
system->total_cap, comm, msg );
|
||||
if( ret != SUCCESS ) {
|
||||
if (ret != SUCCESS) {
|
||||
fprintf( stderr, "no space for workspace: local_cap=%d total_cap=%d",
|
||||
system->local_cap, system->total_cap );
|
||||
fprintf( stderr, "terminating...\n" );
|
||||
|
@ -454,11 +454,11 @@ void ReAllocate( reax_system *system, control_params *control,
|
|||
|
||||
renbr = (data->step - data->prev_steps) % control->reneighbor == 0;
|
||||
/* far neighbors */
|
||||
if( renbr ) {
|
||||
if (renbr) {
|
||||
far_nbrs = *lists + FAR_NBRS;
|
||||
|
||||
if( Nflag || realloc->num_far >= far_nbrs->num_intrs * DANGER_ZONE ) {
|
||||
if( realloc->num_far > far_nbrs->num_intrs ) {
|
||||
if (Nflag || realloc->num_far >= far_nbrs->num_intrs * DANGER_ZONE) {
|
||||
if (realloc->num_far > far_nbrs->num_intrs) {
|
||||
fprintf( stderr, "step%d-ran out of space on far_nbrs: top=%d, max=%d",
|
||||
data->step, realloc->num_far, far_nbrs->num_intrs );
|
||||
MPI_Abort( comm, INSUFFICIENT_MEMORY );
|
||||
|
@ -473,7 +473,7 @@ void ReAllocate( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* hydrogen bonds list */
|
||||
if( control->hbond_cut > 0 ) {
|
||||
if (control->hbond_cut > 0) {
|
||||
Hflag = 0;
|
||||
if( system->numH >= DANGER_ZONE * system->Hcap ||
|
||||
(0 && system->numH <= LOOSE_ZONE * system->Hcap) ) {
|
||||
|
@ -481,7 +481,7 @@ void ReAllocate( reax_system *system, control_params *control,
|
|||
system->Hcap = int(MAX( system->numH * saferzone, mincap ));
|
||||
}
|
||||
|
||||
if( Hflag || realloc->hbonds ) {
|
||||
if (Hflag || realloc->hbonds) {
|
||||
ret = Reallocate_HBonds_List( system, (*lists)+HBONDS, comm );
|
||||
realloc->hbonds = 0;
|
||||
}
|
||||
|
@ -497,10 +497,10 @@ void ReAllocate( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* 3-body list */
|
||||
if( realloc->num_3body > 0 ) {
|
||||
if (realloc->num_3body > 0) {
|
||||
Delete_List( (*lists)+THREE_BODIES, comm );
|
||||
|
||||
if( num_bonds == -1 )
|
||||
if (num_bonds == -1)
|
||||
num_bonds = ((*lists)+BONDS)->num_intrs;
|
||||
|
||||
realloc->num_3body = (int)(MAX(realloc->num_3body*safezone, MIN_3BODIES));
|
||||
|
|
|
@ -110,7 +110,7 @@ void Add_dBond_to_Forces_NPT( int i, int pj, simulation_data *data,
|
|||
/* force */
|
||||
rvec_Add( workspace->f[k], temp );
|
||||
/* pressure */
|
||||
if( k != i ) {
|
||||
if (k != i) {
|
||||
ivec_Sum( rel_box, nbr_k->rel_box, nbr_j->rel_box ); //rel_box(k, i)
|
||||
rvec_iMultiply( ext_press, rel_box, temp );
|
||||
rvec_Add( data->my_ext_press, ext_press );
|
||||
|
@ -189,7 +189,7 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
|
|||
rvec_ScaledAdd( temp, coef.C3dbopi2, workspace->dDeltap_self[i] );
|
||||
rvec_Add( workspace->f[i], temp );
|
||||
|
||||
if( system->pair_ptr->vflag_atom) {
|
||||
if (system->pair_ptr->vflag_atom) {
|
||||
rvec_Scale(fi_tmp, -1.0, temp);
|
||||
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,-1., system->my_atoms[j].x );
|
||||
system->pair_ptr->v_tally(i,fi_tmp,delij);
|
||||
|
@ -208,7 +208,7 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
|
|||
rvec_ScaledAdd( temp, coef.C4dbopi2, workspace->dDeltap_self[j]);
|
||||
rvec_Add( workspace->f[j], temp );
|
||||
|
||||
if( system->pair_ptr->vflag_atom) {
|
||||
if (system->pair_ptr->vflag_atom) {
|
||||
rvec_Scale(fj_tmp, -1.0, temp);
|
||||
rvec_ScaledSum( delji, 1., system->my_atoms[j].x,-1., system->my_atoms[i].x );
|
||||
system->pair_ptr->v_tally(j,fj_tmp,delji);
|
||||
|
@ -225,7 +225,7 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
|
|||
rvec_ScaledAdd( temp, -coef.C3dbopi2, nbr_k->bo_data.dBOp);
|
||||
rvec_Add( workspace->f[k], temp );
|
||||
|
||||
if( system->pair_ptr->vflag_atom ) {
|
||||
if (system->pair_ptr->vflag_atom) {
|
||||
rvec_Scale(fk_tmp, -1.0, temp);
|
||||
rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x);
|
||||
system->pair_ptr->v_tally(k,fk_tmp,delki);
|
||||
|
@ -245,7 +245,7 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
|
|||
rvec_ScaledAdd( temp, -coef.C4dbopi2, nbr_k->bo_data.dBOp);
|
||||
rvec_Add( workspace->f[k], temp );
|
||||
|
||||
if( system->pair_ptr->vflag_atom ) {
|
||||
if (system->pair_ptr->vflag_atom) {
|
||||
rvec_Scale(fk_tmp, -1.0, temp);
|
||||
rvec_ScaledSum(delki,1.,system->my_atoms[k].x,-1.,system->my_atoms[i].x);
|
||||
system->pair_ptr->v_tally(k,fk_tmp,delki);
|
||||
|
@ -271,19 +271,19 @@ int BOp( storage *workspace, reax_list *bonds, double bo_cut,
|
|||
j = nbr_pj->nbr;
|
||||
r2 = SQR(nbr_pj->d);
|
||||
|
||||
if( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0 ) {
|
||||
if (sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0) {
|
||||
C12 = twbp->p_bo1 * pow( nbr_pj->d / twbp->r_s, twbp->p_bo2 );
|
||||
BO_s = (1.0 + bo_cut) * exp( C12 );
|
||||
}
|
||||
else BO_s = C12 = 0.0;
|
||||
|
||||
if( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 ) {
|
||||
if (sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0) {
|
||||
C34 = twbp->p_bo3 * pow( nbr_pj->d / twbp->r_p, twbp->p_bo4 );
|
||||
BO_pi = exp( C34 );
|
||||
}
|
||||
else BO_pi = C34 = 0.0;
|
||||
|
||||
if( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 ) {
|
||||
if (sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0) {
|
||||
C56 = twbp->p_bo5 * pow( nbr_pj->d / twbp->r_pp, twbp->p_bo6 );
|
||||
BO_pi2= exp( C56 );
|
||||
}
|
||||
|
@ -292,7 +292,7 @@ int BOp( storage *workspace, reax_list *bonds, double bo_cut,
|
|||
/* Initially BO values are the uncorrected ones, page 1 */
|
||||
BO = BO_s + BO_pi + BO_pi2;
|
||||
|
||||
if( BO >= bo_cut ) {
|
||||
if (BO >= bo_cut) {
|
||||
/****** bonds i-j and j-i ******/
|
||||
ibond = &( bonds->select.bond_list[btop_i] );
|
||||
btop_j = End_Index( j, bonds );
|
||||
|
@ -410,10 +410,10 @@ void BO( reax_system *system, control_params * /*control*/, simulation_data * /*
|
|||
bo_ij = &( bonds->select.bond_list[pj].bo_data );
|
||||
// fprintf( stderr, "\tj:%d - ubo: %8.3f\n", j+1, bo_ij->BO );
|
||||
|
||||
if( i < j || workspace->bond_mark[j] > 3 ) {
|
||||
if (i < j || workspace->bond_mark[j] > 3) {
|
||||
twbp = &( system->reax_param.tbp[type_i][type_j] );
|
||||
|
||||
if( twbp->ovc < 0.001 && twbp->v13cor < 0.001 ) {
|
||||
if (twbp->ovc < 0.001 && twbp->v13cor < 0.001) {
|
||||
bo_ij->C1dbo = 1.000000;
|
||||
bo_ij->C2dbo = 0.000000;
|
||||
bo_ij->C3dbo = 0.000000;
|
||||
|
@ -435,7 +435,7 @@ void BO( reax_system *system, control_params * /*control*/, simulation_data * /*
|
|||
Deltap_boc_j = workspace->Deltap_boc[j];
|
||||
|
||||
/* on page 1 */
|
||||
if( twbp->ovc >= 0.001 ) {
|
||||
if (twbp->ovc >= 0.001) {
|
||||
/* Correction for overcoordination */
|
||||
exp_p1i = exp( -p_boc1 * Deltap_i );
|
||||
exp_p2i = exp( -p_boc2 * Deltap_i );
|
||||
|
@ -475,7 +475,7 @@ void BO( reax_system *system, control_params * /*control*/, simulation_data * /*
|
|||
Cf1_ij = Cf1_ji = 0.0;
|
||||
}
|
||||
|
||||
if( twbp->v13cor >= 0.001 ) {
|
||||
if (twbp->v13cor >= 0.001) {
|
||||
/* Correction for 1-3 bond orders */
|
||||
exp_f4 =exp(-(twbp->p_boc4 * SQR( bo_ij->BO ) -
|
||||
Deltap_boc_i) * twbp->p_boc3 + twbp->p_boc5);
|
||||
|
@ -527,13 +527,13 @@ void BO( reax_system *system, control_params * /*control*/, simulation_data * /*
|
|||
}
|
||||
|
||||
/* neglect bonds that are < 1e-10 */
|
||||
if( bo_ij->BO < 1e-10 )
|
||||
if (bo_ij->BO < 1e-10)
|
||||
bo_ij->BO = 0.0;
|
||||
if( bo_ij->BO_s < 1e-10 )
|
||||
if (bo_ij->BO_s < 1e-10)
|
||||
bo_ij->BO_s = 0.0;
|
||||
if( bo_ij->BO_pi < 1e-10 )
|
||||
if (bo_ij->BO_pi < 1e-10)
|
||||
bo_ij->BO_pi = 0.0;
|
||||
if( bo_ij->BO_pi2 < 1e-10 )
|
||||
if (bo_ij->BO_pi2 < 1e-10)
|
||||
bo_ij->BO_pi2 = 0.0;
|
||||
|
||||
workspace->total_bond_order[i] += bo_ij->BO; //now keeps total_BO
|
||||
|
@ -576,7 +576,7 @@ void BO( reax_system *system, control_params * /*control*/, simulation_data * /*
|
|||
workspace->Clp[j] = 2.0 * p_lp1 * explp1 * (2.0 + workspace->vlpex[j]);
|
||||
workspace->dDelta_lp[j] = workspace->Clp[j];
|
||||
|
||||
if( sbp_j->mass > 21.0 ) {
|
||||
if (sbp_j->mass > 21.0) {
|
||||
workspace->nlp_temp[j] = 0.5 * (sbp_j->valency_e - sbp_j->valency);
|
||||
workspace->Delta_lp_temp[j] = sbp_j->nlp_opt - workspace->nlp_temp[j];
|
||||
workspace->dDelta_lp_temp[j] = 0.;
|
||||
|
|
|
@ -62,9 +62,9 @@ void Bonds( reax_system *system, control_params * /*control*/,
|
|||
for( pj = start_i; pj < end_i; ++pj ) {
|
||||
j = bonds->select.bond_list[pj].nbr;
|
||||
|
||||
if( system->my_atoms[i].orig_id > system->my_atoms[j].orig_id )
|
||||
if (system->my_atoms[i].orig_id > system->my_atoms[j].orig_id)
|
||||
continue;
|
||||
if( system->my_atoms[i].orig_id == system->my_atoms[j].orig_id ) {
|
||||
if (system->my_atoms[i].orig_id == system->my_atoms[j].orig_id) {
|
||||
if (system->my_atoms[j].x[2] < system->my_atoms[i].x[2]) continue;
|
||||
if (system->my_atoms[j].x[2] == system->my_atoms[i].x[2] &&
|
||||
system->my_atoms[j].x[1] < system->my_atoms[i].x[1]) continue;
|
||||
|
@ -95,7 +95,7 @@ void Bonds( reax_system *system, control_params * /*control*/,
|
|||
-twbp->De_pp * bo_ij->BO_pi2;
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag)
|
||||
if (system->pair_ptr->evflag)
|
||||
system->pair_ptr->ev_tally(i,j,natoms,1,ebond,0.0,0.0,0.0,0.0,0.0);
|
||||
|
||||
/* calculate derivatives of Bond Orders */
|
||||
|
@ -104,7 +104,7 @@ void Bonds( reax_system *system, control_params * /*control*/,
|
|||
bo_ij->Cdbopi2 -= (CEbo + twbp->De_pp);
|
||||
|
||||
/* Stabilisation terminal triple bond */
|
||||
if( bo_ij->BO >= 1.00 ) {
|
||||
if (bo_ij->BO >= 1.00) {
|
||||
if( gp37 == 2 ||
|
||||
(sbp_i->mass == 12.0000 && sbp_j->mass == 15.9990) ||
|
||||
(sbp_j->mass == 12.0000 && sbp_i->mass == 15.9990) ) {
|
||||
|
@ -125,7 +125,7 @@ void Bonds( reax_system *system, control_params * /*control*/,
|
|||
(gp3*exphub1 + 25.0*gp4*exphuov*hulpov*(exphua1+exphub1));
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag)
|
||||
if (system->pair_ptr->evflag)
|
||||
system->pair_ptr->ev_tally(i,j,natoms,1,estriph,0.0,0.0,0.0,0.0,0.0);
|
||||
|
||||
bo_ij->Cdbo += decobdbo;
|
||||
|
|
|
@ -118,16 +118,16 @@ char Read_Control_File( char *control_file, control_params* control,
|
|||
fgets( s, MAX_LINE, fp );
|
||||
Tokenize( s, &tmp );
|
||||
|
||||
if( strcmp(tmp[0], "simulation_name") == 0 ) {
|
||||
if (strcmp(tmp[0], "simulation_name") == 0) {
|
||||
strcpy( control->sim_name, tmp[1] );
|
||||
}
|
||||
else if( strcmp(tmp[0], "ensemble_type") == 0 ) {
|
||||
else if (strcmp(tmp[0], "ensemble_type") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->ensemble = ival;
|
||||
if( ival == iNPT || ival ==sNPT || ival == NPT )
|
||||
if (ival == iNPT || ival ==sNPT || ival == NPT)
|
||||
control->virial = 1;
|
||||
}
|
||||
else if( strcmp(tmp[0], "nsteps") == 0 ) {
|
||||
else if (strcmp(tmp[0], "nsteps") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->nsteps = ival;
|
||||
}
|
||||
|
@ -135,7 +135,7 @@ char Read_Control_File( char *control_file, control_params* control,
|
|||
val = atof(tmp[1]);
|
||||
control->dt = val * 1.e-3; // convert dt from fs to ps!
|
||||
}
|
||||
else if( strcmp(tmp[0], "proc_by_dim") == 0 ) {
|
||||
else if (strcmp(tmp[0], "proc_by_dim") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->procs_by_dim[0] = ival;
|
||||
ival = atoi(tmp[2]);
|
||||
|
@ -146,220 +146,220 @@ char Read_Control_File( char *control_file, control_params* control,
|
|||
control->nprocs = control->procs_by_dim[0]*control->procs_by_dim[1]*
|
||||
control->procs_by_dim[2];
|
||||
}
|
||||
else if( strcmp(tmp[0], "random_vel") == 0 ) {
|
||||
else if (strcmp(tmp[0], "random_vel") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->random_vel = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "restart_format") == 0 ) {
|
||||
else if (strcmp(tmp[0], "restart_format") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->restart_format = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "restart_freq") == 0 ) {
|
||||
else if (strcmp(tmp[0], "restart_freq") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->restart_freq = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "reposition_atoms") == 0 ) {
|
||||
else if (strcmp(tmp[0], "reposition_atoms") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->reposition_atoms = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "restrict_bonds") == 0 ) {
|
||||
else if (strcmp(tmp[0], "restrict_bonds") == 0) {
|
||||
ival = atoi( tmp[1] );
|
||||
control->restrict_bonds = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "remove_CoM_vel") == 0 ) {
|
||||
else if (strcmp(tmp[0], "remove_CoM_vel") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->remove_CoM_vel = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "debug_level") == 0 ) {
|
||||
else if (strcmp(tmp[0], "debug_level") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->debug_level = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "energy_update_freq") == 0 ) {
|
||||
else if (strcmp(tmp[0], "energy_update_freq") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->energy_update_freq = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "reneighbor") == 0 ) {
|
||||
else if (strcmp(tmp[0], "reneighbor") == 0) {
|
||||
ival = atoi( tmp[1] );
|
||||
control->reneighbor = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "vlist_buffer") == 0 ) {
|
||||
else if (strcmp(tmp[0], "vlist_buffer") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->vlist_cut= val + control->nonb_cut;
|
||||
}
|
||||
else if( strcmp(tmp[0], "nbrhood_cutoff") == 0 ) {
|
||||
else if (strcmp(tmp[0], "nbrhood_cutoff") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->bond_cut = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "bond_graph_cutoff") == 0 ) {
|
||||
else if (strcmp(tmp[0], "bond_graph_cutoff") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->bg_cut = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "thb_cutoff") == 0 ) {
|
||||
else if (strcmp(tmp[0], "thb_cutoff") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->thb_cut = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "thb_cutoff_sq") == 0 ) {
|
||||
else if (strcmp(tmp[0], "thb_cutoff_sq") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->thb_cutsq = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "hbond_cutoff") == 0 ) {
|
||||
else if (strcmp(tmp[0], "hbond_cutoff") == 0) {
|
||||
val = atof( tmp[1] );
|
||||
control->hbond_cut = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "ghost_cutoff") == 0 ) {
|
||||
else if (strcmp(tmp[0], "ghost_cutoff") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->user_ghost_cut = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "tabulate_long_range") == 0 ) {
|
||||
else if (strcmp(tmp[0], "tabulate_long_range") == 0) {
|
||||
ival = atoi( tmp[1] );
|
||||
control->tabulate = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "qeq_freq") == 0 ) {
|
||||
else if (strcmp(tmp[0], "qeq_freq") == 0) {
|
||||
ival = atoi( tmp[1] );
|
||||
control->qeq_freq = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "q_err") == 0 ) {
|
||||
else if (strcmp(tmp[0], "q_err") == 0) {
|
||||
val = atof( tmp[1] );
|
||||
control->q_err = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "ilu_refactor") == 0 ) {
|
||||
else if (strcmp(tmp[0], "ilu_refactor") == 0) {
|
||||
ival = atoi( tmp[1] );
|
||||
control->refactor = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "ilu_droptol") == 0 ) {
|
||||
else if (strcmp(tmp[0], "ilu_droptol") == 0) {
|
||||
val = atof( tmp[1] );
|
||||
control->droptol = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "temp_init") == 0 ) {
|
||||
else if (strcmp(tmp[0], "temp_init") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->T_init = val;
|
||||
|
||||
if( control->T_init < 0.1 )
|
||||
if (control->T_init < 0.1)
|
||||
control->T_init = 0.1;
|
||||
}
|
||||
else if( strcmp(tmp[0], "temp_final") == 0 ) {
|
||||
else if (strcmp(tmp[0], "temp_final") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->T_final = val;
|
||||
|
||||
if( control->T_final < 0.1 )
|
||||
if (control->T_final < 0.1)
|
||||
control->T_final = 0.1;
|
||||
}
|
||||
else if( strcmp(tmp[0], "t_mass") == 0 ) {
|
||||
else if (strcmp(tmp[0], "t_mass") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->Tau_T = val * 1.e-3; // convert t_mass from fs to ps
|
||||
}
|
||||
else if( strcmp(tmp[0], "t_mode") == 0 ) {
|
||||
else if (strcmp(tmp[0], "t_mode") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->T_mode = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "t_rate") == 0 ) {
|
||||
else if (strcmp(tmp[0], "t_rate") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->T_rate = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "t_freq") == 0 ) {
|
||||
else if (strcmp(tmp[0], "t_freq") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->T_freq = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "pressure") == 0 ) {
|
||||
if( control->ensemble == iNPT ) {
|
||||
else if (strcmp(tmp[0], "pressure") == 0) {
|
||||
if (control->ensemble == iNPT) {
|
||||
control->P[0] = control->P[1] = control->P[2] = atof(tmp[1]);
|
||||
}
|
||||
else if( control->ensemble == sNPT ) {
|
||||
else if (control->ensemble == sNPT) {
|
||||
control->P[0] = atof(tmp[1]);
|
||||
control->P[1] = atof(tmp[2]);
|
||||
control->P[2] = atof(tmp[3]);
|
||||
}
|
||||
}
|
||||
else if( strcmp(tmp[0], "p_mass") == 0 ) {
|
||||
else if (strcmp(tmp[0], "p_mass") == 0) {
|
||||
// convert p_mass from fs to ps
|
||||
if( control->ensemble == iNPT ) {
|
||||
if (control->ensemble == iNPT) {
|
||||
control->Tau_P[0] = control->Tau_P[1] = control->Tau_P[2] =
|
||||
atof(tmp[1]) * 1.e-3;
|
||||
}
|
||||
else if( control->ensemble == sNPT ) {
|
||||
else if (control->ensemble == sNPT) {
|
||||
control->Tau_P[0] = atof(tmp[1]) * 1.e-3;
|
||||
control->Tau_P[1] = atof(tmp[2]) * 1.e-3;
|
||||
control->Tau_P[2] = atof(tmp[3]) * 1.e-3;
|
||||
}
|
||||
}
|
||||
else if( strcmp(tmp[0], "pt_mass") == 0 ) {
|
||||
else if (strcmp(tmp[0], "pt_mass") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->Tau_PT[0] = control->Tau_PT[1] = control->Tau_PT[2] =
|
||||
val * 1.e-3; // convert pt_mass from fs to ps
|
||||
}
|
||||
else if( strcmp(tmp[0], "compress") == 0 ) {
|
||||
else if (strcmp(tmp[0], "compress") == 0) {
|
||||
val = atof(tmp[1]);
|
||||
control->compressibility = val;
|
||||
}
|
||||
else if( strcmp(tmp[0], "press_mode") == 0 ) {
|
||||
else if (strcmp(tmp[0], "press_mode") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->press_mode = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "geo_format") == 0 ) {
|
||||
else if (strcmp(tmp[0], "geo_format") == 0) {
|
||||
ival = atoi( tmp[1] );
|
||||
control->geo_format = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "write_freq") == 0 ) {
|
||||
else if (strcmp(tmp[0], "write_freq") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->write_steps = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "traj_compress") == 0 ) {
|
||||
else if (strcmp(tmp[0], "traj_compress") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->traj_compress = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "traj_method") == 0 ) {
|
||||
else if (strcmp(tmp[0], "traj_method") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->traj_method = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "traj_title") == 0 ) {
|
||||
else if (strcmp(tmp[0], "traj_title") == 0) {
|
||||
strcpy( out_control->traj_title, tmp[1] );
|
||||
}
|
||||
else if( strcmp(tmp[0], "atom_info") == 0 ) {
|
||||
else if (strcmp(tmp[0], "atom_info") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->atom_info += ival * 4;
|
||||
}
|
||||
else if( strcmp(tmp[0], "atom_velocities") == 0 ) {
|
||||
else if (strcmp(tmp[0], "atom_velocities") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->atom_info += ival * 2;
|
||||
}
|
||||
else if( strcmp(tmp[0], "atom_forces") == 0 ) {
|
||||
else if (strcmp(tmp[0], "atom_forces") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->atom_info += ival * 1;
|
||||
}
|
||||
else if( strcmp(tmp[0], "bond_info") == 0 ) {
|
||||
else if (strcmp(tmp[0], "bond_info") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->bond_info = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "angle_info") == 0 ) {
|
||||
else if (strcmp(tmp[0], "angle_info") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
out_control->angle_info = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "molecular_analysis") == 0 ) {
|
||||
else if (strcmp(tmp[0], "molecular_analysis") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->molecular_analysis = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "ignore") == 0 ) {
|
||||
else if (strcmp(tmp[0], "ignore") == 0) {
|
||||
control->num_ignored = atoi(tmp[1]);
|
||||
for( i = 0; i < control->num_ignored; ++i )
|
||||
control->ignore[atoi(tmp[i+2])] = 1;
|
||||
}
|
||||
else if( strcmp(tmp[0], "dipole_anal") == 0 ) {
|
||||
else if (strcmp(tmp[0], "dipole_anal") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->dipole_anal = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "freq_dipole_anal") == 0 ) {
|
||||
else if (strcmp(tmp[0], "freq_dipole_anal") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->freq_dipole_anal = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "diffusion_coef") == 0 ) {
|
||||
else if (strcmp(tmp[0], "diffusion_coef") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->diffusion_coef = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "freq_diffusion_coef") == 0 ) {
|
||||
else if (strcmp(tmp[0], "freq_diffusion_coef") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->freq_diffusion_coef = ival;
|
||||
}
|
||||
else if( strcmp(tmp[0], "restrict_type") == 0 ) {
|
||||
else if (strcmp(tmp[0], "restrict_type") == 0) {
|
||||
ival = atoi(tmp[1]);
|
||||
control->restrict_type = ival;
|
||||
}
|
||||
|
@ -370,7 +370,7 @@ char Read_Control_File( char *control_file, control_params* control,
|
|||
}
|
||||
|
||||
/* determine target T */
|
||||
if( control->T_mode == 0 )
|
||||
if (control->T_mode == 0)
|
||||
control->T = control->T_final;
|
||||
else control->T = control->T_init;
|
||||
|
||||
|
|
|
@ -260,7 +260,7 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
|
|||
|
||||
if (reax->sbp[i].rcore2>0.01 && reax->sbp[i].acore2>0.01) { // Inner-wall
|
||||
if (reax->sbp[i].gamma_w>0.5) { // Shielding vdWaals
|
||||
if( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 3 ) {
|
||||
if (reax->gp.vdw_type != 0 && reax->gp.vdw_type != 3) {
|
||||
if (errorflag && (me == 0))
|
||||
fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n" \
|
||||
"Force field parameters for element %s\n" \
|
||||
|
@ -274,7 +274,7 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
|
|||
reax->gp.vdw_type = 3;
|
||||
}
|
||||
} else { // No shielding vdWaals parameters present
|
||||
if( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 2 ) {
|
||||
if (reax->gp.vdw_type != 0 && reax->gp.vdw_type != 2) {
|
||||
if (me == 0)
|
||||
fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n" \
|
||||
"Force field parameters for element %s\n" \
|
||||
|
@ -289,7 +289,7 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
|
|||
}
|
||||
} else { // No Inner wall parameters present
|
||||
if (reax->sbp[i].gamma_w>0.5) { // Shielding vdWaals
|
||||
if( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 1 ) {
|
||||
if (reax->gp.vdw_type != 0 && reax->gp.vdw_type != 1) {
|
||||
if (me == 0)
|
||||
fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n" \
|
||||
"Force field parameters for element %s\n" \
|
||||
|
@ -641,7 +641,7 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
|
|||
reax->fbp[n][m][k][j].prm[0].p_cot1 = val;
|
||||
}
|
||||
} else { /* This means the entry is of the form 0-X-Y-0 */
|
||||
if( k < reax->num_atom_types && m < reax->num_atom_types )
|
||||
if (k < reax->num_atom_types && m < reax->num_atom_types)
|
||||
for( p = 0; p < reax->num_atom_types; p++ )
|
||||
for( o = 0; o < reax->num_atom_types; o++ ) {
|
||||
reax->fbp[p][k][m][o].cnt = 1;
|
||||
|
@ -687,7 +687,7 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
|
|||
m = atoi(tmp[2]) - 1;
|
||||
|
||||
|
||||
if( j < reax->num_atom_types && m < reax->num_atom_types ) {
|
||||
if (j < reax->num_atom_types && m < reax->num_atom_types) {
|
||||
val = atof(tmp[3]);
|
||||
reax->hbp[j][k][m].r0_hb = val;
|
||||
|
||||
|
|
|
@ -55,7 +55,7 @@ void Init_Force_Functions( control_params *control )
|
|||
Interaction_Functions[2] = Atom_Energy; //Dummy_Interaction;
|
||||
Interaction_Functions[3] = Valence_Angles; //Dummy_Interaction;
|
||||
Interaction_Functions[4] = Torsion_Angles; //Dummy_Interaction;
|
||||
if( control->hbond_cut > 0 )
|
||||
if (control->hbond_cut > 0)
|
||||
Interaction_Functions[5] = Hydrogen_Bonds;
|
||||
else Interaction_Functions[5] = Dummy_Interaction;
|
||||
Interaction_Functions[6] = Dummy_Interaction; //empty
|
||||
|
@ -87,7 +87,7 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
|
|||
{
|
||||
|
||||
/* van der Waals and Coulomb interactions */
|
||||
if( control->tabulate == 0 )
|
||||
if (control->tabulate == 0)
|
||||
vdW_Coulomb_Energy( system, control, data, workspace,
|
||||
lists, out_control );
|
||||
else
|
||||
|
@ -105,8 +105,8 @@ void Compute_Total_Force( reax_system *system, control_params *control,
|
|||
|
||||
for( i = 0; i < system->N; ++i )
|
||||
for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
|
||||
if( i < bonds->select.bond_list[pj].nbr ) {
|
||||
if( control->virial == 0 )
|
||||
if (i < bonds->select.bond_list[pj].nbr) {
|
||||
if (control->virial == 0)
|
||||
Add_dBond_to_Forces( system, i, pj, workspace, lists );
|
||||
else
|
||||
Add_dBond_to_Forces_NPT( i, pj, data, workspace, lists );
|
||||
|
@ -123,17 +123,17 @@ void Validate_Lists( reax_system *system, storage * /*workspace*/, reax_list **l
|
|||
double saferzone = system->saferzone;
|
||||
|
||||
/* bond list */
|
||||
if( N > 0 ) {
|
||||
if (N > 0) {
|
||||
bonds = *lists + BONDS;
|
||||
|
||||
for( i = 0; i < N; ++i ) {
|
||||
system->my_atoms[i].num_bonds = MAX(Num_Entries(i,bonds)*2, MIN_BONDS);
|
||||
|
||||
if( i < N-1 )
|
||||
if (i < N-1)
|
||||
comp = Start_Index(i+1, bonds);
|
||||
else comp = bonds->num_intrs;
|
||||
|
||||
if( End_Index(i, bonds) > comp ) {
|
||||
if (End_Index(i, bonds) > comp) {
|
||||
fprintf( stderr, "step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
|
||||
step, i, End_Index(i,bonds), comp );
|
||||
MPI_Abort( comm, INSUFFICIENT_MEMORY );
|
||||
|
@ -143,12 +143,12 @@ void Validate_Lists( reax_system *system, storage * /*workspace*/, reax_list **l
|
|||
|
||||
|
||||
/* hbonds list */
|
||||
if( numH > 0 ) {
|
||||
if (numH > 0) {
|
||||
hbonds = *lists + HBONDS;
|
||||
|
||||
for( i = 0; i < N; ++i ) {
|
||||
Hindex = system->my_atoms[i].Hindex;
|
||||
if( Hindex > -1 ) {
|
||||
if (Hindex > -1) {
|
||||
system->my_atoms[i].num_hbonds =
|
||||
(int)(MAX( Num_Entries(Hindex, hbonds)*saferzone, MIN_HBONDS ));
|
||||
|
||||
|
@ -156,11 +156,11 @@ void Validate_Lists( reax_system *system, storage * /*workspace*/, reax_list **l
|
|||
//(Start_Index(i+1,hbonds)-Start_Index(i,hbonds))*0.90/*DANGER_ZONE*/){
|
||||
// workspace->realloc.hbonds = 1;
|
||||
|
||||
if( Hindex < numH-1 )
|
||||
if (Hindex < numH-1)
|
||||
comp = Start_Index(Hindex+1, hbonds);
|
||||
else comp = hbonds->num_intrs;
|
||||
|
||||
if( End_Index(Hindex, hbonds) > comp ) {
|
||||
if (End_Index(Hindex, hbonds) > comp) {
|
||||
fprintf(stderr,"step%d-hbondchk failed: H=%d end(H)=%d str(H+1)=%d\n",
|
||||
step, Hindex, End_Index(Hindex,hbonds), comp );
|
||||
MPI_Abort( comm, INSUFFICIENT_MEMORY );
|
||||
|
@ -212,7 +212,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
|
|||
btop_i = End_Index( i, bonds );
|
||||
sbp_i = &(system->reax_param.sbp[type_i]);
|
||||
|
||||
if( i < system->n ) {
|
||||
if (i < system->n) {
|
||||
local = 1;
|
||||
cutoff = MAX( control->hbond_cut, control->bond_cut );
|
||||
} else {
|
||||
|
@ -222,9 +222,9 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
|
|||
|
||||
ihb = -1;
|
||||
ihb_top = -1;
|
||||
if( local && control->hbond_cut > 0 ) {
|
||||
if (local && control->hbond_cut > 0) {
|
||||
ihb = sbp_i->p_hbond;
|
||||
if( ihb == 1 )
|
||||
if (ihb == 1)
|
||||
ihb_top = End_Index( atom_i->Hindex, hbonds );
|
||||
else ihb_top = -1;
|
||||
}
|
||||
|
@ -235,8 +235,8 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
|
|||
j = nbr_pj->nbr;
|
||||
atom_j = &(system->my_atoms[j]);
|
||||
|
||||
if( renbr ) {
|
||||
if( nbr_pj->d <= cutoff )
|
||||
if (renbr) {
|
||||
if (nbr_pj->d <= cutoff)
|
||||
flag = 1;
|
||||
else flag = 0;
|
||||
} else {
|
||||
|
@ -244,7 +244,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
|
|||
nbr_pj->dvec[1] = atom_j->x[1] - atom_i->x[1];
|
||||
nbr_pj->dvec[2] = atom_j->x[2] - atom_i->x[2];
|
||||
nbr_pj->d = rvec_Norm_Sqr( nbr_pj->dvec );
|
||||
if( nbr_pj->d <= SQR(cutoff) ) {
|
||||
if (nbr_pj->d <= SQR(cutoff)) {
|
||||
nbr_pj->d = sqrt(nbr_pj->d);
|
||||
flag = 1;
|
||||
} else {
|
||||
|
@ -252,26 +252,26 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
|
|||
}
|
||||
}
|
||||
|
||||
if( flag ) {
|
||||
if (flag) {
|
||||
type_j = atom_j->type;
|
||||
if (type_j < 0) continue;
|
||||
sbp_j = &(system->reax_param.sbp[type_j]);
|
||||
twbp = &(system->reax_param.tbp[type_i][type_j]);
|
||||
|
||||
if( local ) {
|
||||
if (local) {
|
||||
/* hydrogen bond lists */
|
||||
if( control->hbond_cut > 0 && (ihb==1 || ihb==2) &&
|
||||
if (control->hbond_cut > 0 && (ihb==1 || ihb==2) &&
|
||||
nbr_pj->d <= control->hbond_cut ) {
|
||||
// fprintf( stderr, "%d %d\n", atom1, atom2 );
|
||||
jhb = sbp_j->p_hbond;
|
||||
if( ihb == 1 && jhb == 2 ) {
|
||||
if (ihb == 1 && jhb == 2) {
|
||||
hbonds->select.hbond_list[ihb_top].nbr = j;
|
||||
hbonds->select.hbond_list[ihb_top].scl = 1;
|
||||
hbonds->select.hbond_list[ihb_top].ptr = nbr_pj;
|
||||
++ihb_top;
|
||||
++num_hbonds;
|
||||
}
|
||||
else if( j < system->n && ihb == 2 && jhb == 1 ) {
|
||||
else if (j < system->n && ihb == 2 && jhb == 1) {
|
||||
jhb_top = End_Index( atom_j->Hindex, hbonds );
|
||||
hbonds->select.hbond_list[jhb_top].nbr = i;
|
||||
hbonds->select.hbond_list[jhb_top].scl = -1;
|
||||
|
@ -282,16 +282,16 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
|
|||
}
|
||||
}
|
||||
|
||||
if( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
|
||||
if (//(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
|
||||
nbr_pj->d <= control->bond_cut &&
|
||||
BOp( workspace, bonds, control->bo_cut,
|
||||
i , btop_i, nbr_pj, sbp_i, sbp_j, twbp ) ) {
|
||||
num_bonds += 2;
|
||||
++btop_i;
|
||||
|
||||
if( workspace->bond_mark[j] > workspace->bond_mark[i] + 1 )
|
||||
if (workspace->bond_mark[j] > workspace->bond_mark[i] + 1)
|
||||
workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
|
||||
else if( workspace->bond_mark[i] > workspace->bond_mark[j] + 1 ) {
|
||||
else if (workspace->bond_mark[i] > workspace->bond_mark[j] + 1) {
|
||||
workspace->bond_mark[i] = workspace->bond_mark[j] + 1;
|
||||
}
|
||||
}
|
||||
|
@ -299,7 +299,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
Set_End_Index( i, btop_i, bonds );
|
||||
if( local && ihb == 1 )
|
||||
if (local && ihb == 1)
|
||||
Set_End_Index( atom_i->Hindex, ihb_top, hbonds );
|
||||
}
|
||||
|
||||
|
@ -349,7 +349,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
|
|||
end_i = End_Index(i, far_nbrs);
|
||||
sbp_i = &(system->reax_param.sbp[type_i]);
|
||||
|
||||
if( i < system->n ) {
|
||||
if (i < system->n) {
|
||||
local = 1;
|
||||
cutoff = control->nonb_cut;
|
||||
++(*Htop);
|
||||
|
@ -372,15 +372,15 @@ void Estimate_Storages( reax_system *system, control_params *control,
|
|||
sbp_j = &(system->reax_param.sbp[type_j]);
|
||||
twbp = &(system->reax_param.tbp[type_i][type_j]);
|
||||
|
||||
if( local ) {
|
||||
if( j < system->n || atom_i->orig_id < atom_j->orig_id ) //tryQEq ||1
|
||||
if (local) {
|
||||
if (j < system->n || atom_i->orig_id < atom_j->orig_id) //tryQEq ||1
|
||||
++(*Htop);
|
||||
|
||||
/* hydrogen bond lists */
|
||||
if( control->hbond_cut > 0.1 && (ihb==1 || ihb==2) &&
|
||||
if (control->hbond_cut > 0.1 && (ihb==1 || ihb==2) &&
|
||||
nbr_pj->d <= control->hbond_cut ) {
|
||||
jhb = sbp_j->p_hbond;
|
||||
if( ihb == 1 && jhb == 2 )
|
||||
if (ihb == 1 && jhb == 2)
|
||||
++hb_top[i];
|
||||
else if( j < system->n && ihb == 2 && jhb == 1 )
|
||||
++hb_top[j];
|
||||
|
@ -388,20 +388,20 @@ void Estimate_Storages( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* uncorrected bond orders */
|
||||
if( nbr_pj->d <= control->bond_cut ) {
|
||||
if( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0) {
|
||||
if (nbr_pj->d <= control->bond_cut) {
|
||||
if (sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0) {
|
||||
C12 = twbp->p_bo1 * pow( r_ij / twbp->r_s, twbp->p_bo2 );
|
||||
BO_s = (1.0 + control->bo_cut) * exp( C12 );
|
||||
}
|
||||
else BO_s = C12 = 0.0;
|
||||
|
||||
if( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0) {
|
||||
if (sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0) {
|
||||
C34 = twbp->p_bo3 * pow( r_ij / twbp->r_p, twbp->p_bo4 );
|
||||
BO_pi = exp( C34 );
|
||||
}
|
||||
else BO_pi = C34 = 0.0;
|
||||
|
||||
if( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0) {
|
||||
if (sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0) {
|
||||
C56 = twbp->p_bo5 * pow( r_ij / twbp->r_pp, twbp->p_bo6 );
|
||||
BO_pi2= exp( C56 );
|
||||
}
|
||||
|
@ -410,7 +410,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
|
|||
/* Initially BO values are the uncorrected ones, page 1 */
|
||||
BO = BO_s + BO_pi + BO_pi2;
|
||||
|
||||
if( BO >= control->bo_cut ) {
|
||||
if (BO >= control->bo_cut) {
|
||||
++bond_top[i];
|
||||
++bond_top[j];
|
||||
}
|
||||
|
|
|
@ -63,7 +63,7 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
|
|||
hbond_list = hbonds->select.hbond_list;
|
||||
|
||||
for( j = 0; j < system->n; ++j )
|
||||
if( system->reax_param.sbp[system->my_atoms[j].type].p_hbond == 1 ) {
|
||||
if (system->reax_param.sbp[system->my_atoms[j].type].p_hbond == 1) {
|
||||
type_j = system->my_atoms[j].type;
|
||||
start_j = Start_Index(j, bonds);
|
||||
end_j = End_Index(j, bonds);
|
||||
|
@ -98,7 +98,7 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
|
|||
pbond_ij = &( bonds->select.bond_list[pi] );
|
||||
i = pbond_ij->nbr;
|
||||
|
||||
if( system->my_atoms[i].orig_id != system->my_atoms[k].orig_id ) {
|
||||
if (system->my_atoms[i].orig_id != system->my_atoms[k].orig_id) {
|
||||
bo_ij = &(pbond_ij->bo_data);
|
||||
type_i = system->my_atoms[i].type;
|
||||
if (type_i < 0) continue;
|
||||
|
@ -133,7 +133,7 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
|
|||
/* hydrogen bond forces */
|
||||
bo_ij->Cdbo += CEhb1; // dbo term
|
||||
|
||||
if( control->virial == 0 ) {
|
||||
if (control->virial == 0) {
|
||||
// dcos terms
|
||||
rvec_ScaledAdd( workspace->f[i], +CEhb2, dcos_theta_di );
|
||||
rvec_ScaledAdd( workspace->f[j], +CEhb2, dcos_theta_dj );
|
||||
|
|
|
@ -52,7 +52,7 @@ int Init_System( reax_system *system, control_params *control, char * /*msg*/ )
|
|||
|
||||
/* estimate numH and Hcap */
|
||||
system->numH = 0;
|
||||
if( control->hbond_cut > 0 )
|
||||
if (control->hbond_cut > 0)
|
||||
for( i = 0; i < system->n; ++i ) {
|
||||
atom = &(system->my_atoms[i]);
|
||||
if (system->reax_param.sbp[ atom->type ].p_hbond == 1 && atom->type >= 0)
|
||||
|
@ -71,7 +71,7 @@ int Init_Simulation_Data( reax_system *system, control_params *control,
|
|||
Reset_Simulation_Data( data, control->virial );
|
||||
|
||||
/* initialize the timer(s) */
|
||||
if( system->my_rank == MASTER_NODE ) {
|
||||
if (system->my_rank == MASTER_NODE) {
|
||||
data->timing.start = Get_Time( );
|
||||
}
|
||||
|
||||
|
@ -89,10 +89,10 @@ void Init_Taper( control_params *control, storage *workspace, MPI_Comm comm )
|
|||
swa = control->nonb_low;
|
||||
swb = control->nonb_cut;
|
||||
|
||||
if( fabs( swa ) > 0.01 )
|
||||
if (fabs( swa ) > 0.01)
|
||||
fprintf( stderr, "Warning: non-zero lower Taper-radius cutoff\n" );
|
||||
|
||||
if( swb < 0 ) {
|
||||
if (swb < 0) {
|
||||
fprintf( stderr, "Negative upper Taper-radius cutoff\n" );
|
||||
MPI_Abort( comm, INVALID_INPUT );
|
||||
}
|
||||
|
@ -125,7 +125,7 @@ int Init_Workspace( reax_system *system, control_params *control,
|
|||
|
||||
ret = Allocate_Workspace( system, control, workspace,
|
||||
system->local_cap, system->total_cap, comm, msg );
|
||||
if( ret != SUCCESS )
|
||||
if (ret != SUCCESS)
|
||||
return ret;
|
||||
|
||||
memset( &(workspace->realloc), 0, sizeof(reallocate_data) );
|
||||
|
@ -168,7 +168,7 @@ int Init_Lists( reax_system *system, control_params *control,
|
|||
Estimate_Storages( system, control, lists,
|
||||
&Htop, hb_top, bond_top, &num_3body, comm );
|
||||
|
||||
if( control->hbond_cut > 0 ) {
|
||||
if (control->hbond_cut > 0) {
|
||||
/* init H indexes */
|
||||
total_hbonds = 0;
|
||||
for( i = 0; i < system->n; ++i ) {
|
||||
|
@ -219,7 +219,7 @@ void Initialize( reax_system *system, control_params *control,
|
|||
char msg[MAX_STR];
|
||||
|
||||
|
||||
if( Init_MPI_Datatypes(system, workspace, mpi_data, comm, msg) == FAILURE ) {
|
||||
if (Init_MPI_Datatypes(system, workspace, mpi_data, comm, msg) == FAILURE) {
|
||||
fprintf( stderr, "p%d: init_mpi_datatypes: could not create datatypes\n",
|
||||
system->my_rank );
|
||||
fprintf( stderr, "p%d: mpi_data couldn't be initialized! terminating.\n",
|
||||
|
@ -234,15 +234,15 @@ void Initialize( reax_system *system, control_params *control,
|
|||
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
|
||||
}
|
||||
|
||||
if( Init_Simulation_Data( system, control, data, msg ) == FAILURE ) {
|
||||
if (Init_Simulation_Data( system, control, data, msg ) == FAILURE) {
|
||||
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
|
||||
fprintf( stderr, "p%d: sim_data couldn't be initialized! terminating.\n",
|
||||
system->my_rank );
|
||||
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
|
||||
}
|
||||
|
||||
if( Init_Workspace( system, control, workspace, mpi_data->world, msg ) ==
|
||||
FAILURE ) {
|
||||
if (Init_Workspace( system, control, workspace, mpi_data->world, msg ) ==
|
||||
FAILURE) {
|
||||
fprintf( stderr, "p%d:init_workspace: not enough memory\n",
|
||||
system->my_rank );
|
||||
fprintf( stderr, "p%d:workspace couldn't be initialized! terminating.\n",
|
||||
|
@ -250,23 +250,23 @@ void Initialize( reax_system *system, control_params *control,
|
|||
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
|
||||
}
|
||||
|
||||
if( Init_Lists( system, control, data, workspace, lists, mpi_data, msg ) ==
|
||||
FAILURE ) {
|
||||
if (Init_Lists( system, control, data, workspace, lists, mpi_data, msg ) ==
|
||||
FAILURE) {
|
||||
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
|
||||
fprintf( stderr, "p%d: system could not be initialized! terminating.\n",
|
||||
system->my_rank );
|
||||
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
|
||||
}
|
||||
|
||||
if( Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE) {
|
||||
if (Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE) {
|
||||
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
|
||||
fprintf( stderr, "p%d: could not open output files! terminating...\n",
|
||||
system->my_rank );
|
||||
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
|
||||
}
|
||||
|
||||
if( control->tabulate ) {
|
||||
if( Init_Lookup_Tables( system, control, workspace, mpi_data, msg ) == FAILURE ) {
|
||||
if (control->tabulate) {
|
||||
if (Init_Lookup_Tables( system, control, workspace, mpi_data, msg ) == FAILURE) {
|
||||
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
|
||||
fprintf( stderr, "p%d: couldn't create lookup table! terminating.\n",
|
||||
system->my_rank );
|
||||
|
|
|
@ -43,17 +43,17 @@ int Init_Output_Files( reax_system *system, control_params *control,
|
|||
|
||||
if (out_control->write_steps > 0) {
|
||||
ret = Init_Traj( system, control, out_control, mpi_data, msg );
|
||||
if( ret == FAILURE )
|
||||
if (ret == FAILURE)
|
||||
return ret;
|
||||
}
|
||||
|
||||
if( system->my_rank == MASTER_NODE ) {
|
||||
if (system->my_rank == MASTER_NODE) {
|
||||
/* These files are written only by the master node */
|
||||
if( out_control->energy_update_freq > 0 ) {
|
||||
if (out_control->energy_update_freq > 0) {
|
||||
|
||||
/* init potentials file */
|
||||
sprintf( temp, "%s.pot", control->sim_name );
|
||||
if( (out_control->pot = fopen( temp, "w" )) != NULL ) {
|
||||
if ((out_control->pot = fopen( temp, "w" )) != NULL) {
|
||||
fflush( out_control->pot );
|
||||
}
|
||||
else {
|
||||
|
@ -69,7 +69,7 @@ int Init_Output_Files( reax_system *system, control_params *control,
|
|||
control->ensemble == iNPT ||
|
||||
control->ensemble == sNPT ) {
|
||||
sprintf( temp, "%s.prs", control->sim_name );
|
||||
if( (out_control->prs = fopen( temp, "w" )) != NULL ) {
|
||||
if ((out_control->prs = fopen( temp, "w" )) != NULL) {
|
||||
fprintf(out_control->prs,"%8s%13s%13s%13s%13s%13s%13s%13s\n",
|
||||
"step", "Pint/norm[x]", "Pint/norm[y]", "Pint/norm[z]",
|
||||
"Pext/Ptot[x]", "Pext/Ptot[y]", "Pext/Ptot[z]", "Pkin/V" );
|
||||
|
@ -90,11 +90,11 @@ int Init_Output_Files( reax_system *system, control_params *control,
|
|||
int Close_Output_Files( reax_system *system, control_params *control,
|
||||
output_controls *out_control, mpi_datatypes * /*mpi_data*/ )
|
||||
{
|
||||
if( out_control->write_steps > 0 )
|
||||
if (out_control->write_steps > 0)
|
||||
End_Traj( system->my_rank, out_control );
|
||||
|
||||
if( system->my_rank == MASTER_NODE ) {
|
||||
if( out_control->energy_update_freq > 0 ) {
|
||||
if (system->my_rank == MASTER_NODE) {
|
||||
if (out_control->energy_update_freq > 0) {
|
||||
fclose( out_control->pot );
|
||||
}
|
||||
|
||||
|
|
|
@ -97,7 +97,7 @@ int Make_List(int n, int num_intrs, int type, reax_list *l, MPI_Comm comm)
|
|||
|
||||
void Delete_List( reax_list *l, MPI_Comm comm )
|
||||
{
|
||||
if( l->allocated == 0 )
|
||||
if (l->allocated == 0)
|
||||
return;
|
||||
l->allocated = 0;
|
||||
|
||||
|
|
|
@ -198,9 +198,9 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
|
|||
MPI_INT, MPI_SUM, mpi_data->world );
|
||||
|
||||
for( i = 0; i < num_atom_types; ++i ) {
|
||||
if( aggregated[i] ) {
|
||||
if (aggregated[i]) {
|
||||
for( j = i; j < num_atom_types; ++j ) {
|
||||
if( aggregated[j] ) {
|
||||
if (aggregated[j]) {
|
||||
LR[i][j].xmin = 0;
|
||||
LR[i][j].xmax = control->nonb_cut;
|
||||
LR[i][j].n = control->tabulate + 2;
|
||||
|
@ -290,7 +290,7 @@ void Deallocate_Lookup_Tables( reax_system *system )
|
|||
|
||||
for( i = 0; i < ntypes; ++i ) {
|
||||
for( j = i; j < ntypes; ++j )
|
||||
if( LR[i][j].n ) {
|
||||
if (LR[i][j].n) {
|
||||
sfree( LR[i][j].y, "LR[i,j].y" );
|
||||
sfree( LR[i][j].H, "LR[i,j].H" );
|
||||
sfree( LR[i][j].vdW, "LR[i,j].vdW" );
|
||||
|
|
|
@ -91,23 +91,23 @@ void Atom_Energy( reax_system *system, control_params *control,
|
|||
workspace->CdDelta[i] += CElp; // lp - 1st term
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag)
|
||||
if (system->pair_ptr->evflag)
|
||||
system->pair_ptr->ev_tally(i,i,system->n,1,e_lp,0.0,0.0,0.0,0.0,0.0);
|
||||
|
||||
/* correction for C2 */
|
||||
if( p_lp3 > 0.001 && !strcmp(system->reax_param.sbp[type_i].name, "C") )
|
||||
if (p_lp3 > 0.001 && !strcmp(system->reax_param.sbp[type_i].name, "C"))
|
||||
for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj ) {
|
||||
j = bonds->select.bond_list[pj].nbr;
|
||||
type_j = system->my_atoms[j].type;
|
||||
if (type_j < 0) continue;
|
||||
|
||||
if( !strcmp( system->reax_param.sbp[type_j].name, "C" ) ) {
|
||||
if (!strcmp( system->reax_param.sbp[type_j].name, "C" )) {
|
||||
twbp = &( system->reax_param.tbp[type_i][type_j]);
|
||||
bo_ij = &( bonds->select.bond_list[pj].bo_data );
|
||||
Di = workspace->Delta[i];
|
||||
vov3 = bo_ij->BO - Di - 0.040*pow(Di, 4.);
|
||||
|
||||
if( vov3 > 3. ) {
|
||||
if (vov3 > 3.) {
|
||||
data->my_en.e_lp += e_lph = p_lp3 * SQR(vov3-3.0);
|
||||
|
||||
deahu2dbo = 2.*p_lp3*(vov3 - 3.);
|
||||
|
@ -117,7 +117,7 @@ void Atom_Energy( reax_system *system, control_params *control,
|
|||
workspace->CdDelta[i] += deahu2dsbo;
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag)
|
||||
if (system->pair_ptr->evflag)
|
||||
system->pair_ptr->ev_tally(i,j,system->n,1,e_lph,0.0,0.0,0.0,0.0,0.0);
|
||||
|
||||
}
|
||||
|
@ -132,7 +132,7 @@ void Atom_Energy( reax_system *system, control_params *control,
|
|||
sbp_i = &(system->reax_param.sbp[ type_i ]);
|
||||
|
||||
/* over-coordination energy */
|
||||
if( sbp_i->mass > 21.0 )
|
||||
if (sbp_i->mass > 21.0)
|
||||
dfvl = 0.0;
|
||||
else dfvl = 1.0; // only for 1st-row elements
|
||||
|
||||
|
@ -201,7 +201,7 @@ void Atom_Energy( reax_system *system, control_params *control,
|
|||
p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1) + CEunder2;
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag) {
|
||||
if (system->pair_ptr->evflag) {
|
||||
eng_tmp = e_ov;
|
||||
if (numbonds > 0 || control->enobondsflag)
|
||||
eng_tmp += e_un;
|
||||
|
|
|
@ -172,7 +172,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
|
|||
( dTap - Tap * r_ij / dr3gamij_1 ) / dr3gamij_3;
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
|
||||
if (system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
|
||||
pe_vdw = Tap * (e_vdW + e_core + e_lg);
|
||||
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
|
||||
-1., system->my_atoms[j].x );
|
||||
|
@ -181,7 +181,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
|
|||
f_tmp,delij[0],delij[1],delij[2]);
|
||||
}
|
||||
|
||||
if( control->virial == 0 ) {
|
||||
if (control->virial == 0) {
|
||||
rvec_ScaledAdd( workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec );
|
||||
rvec_ScaledAdd( workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec );
|
||||
} else { /* NPT, iNPT or sNPT */
|
||||
|
@ -263,7 +263,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
|
|||
|
||||
/* Cubic Spline Interpolation */
|
||||
r = (int)(r_ij * t->inv_dx);
|
||||
if( r == 0 ) ++r;
|
||||
if (r == 0) ++r;
|
||||
base = (double)(r+1) * t->dx;
|
||||
dif = r_ij - base;
|
||||
|
||||
|
@ -285,7 +285,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
|
|||
CEclmb *= system->my_atoms[i].q * system->my_atoms[j].q;
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
|
||||
if (system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
|
||||
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
|
||||
-1., system->my_atoms[j].x );
|
||||
f_tmp = -(CEvd + CEclmb);
|
||||
|
@ -293,7 +293,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
|
|||
f_tmp,delij[0],delij[1],delij[2]);
|
||||
}
|
||||
|
||||
if( control->virial == 0 ) {
|
||||
if (control->virial == 0) {
|
||||
rvec_ScaledAdd( workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec );
|
||||
rvec_ScaledAdd( workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec );
|
||||
} else { // NPT, iNPT or sNPT
|
||||
|
@ -330,7 +330,7 @@ void Compute_Polarization_Energy( reax_system *system, simulation_data *data )
|
|||
data->my_en.e_pol += en_tmp;
|
||||
|
||||
/* tally into per-atom energy */
|
||||
if( system->pair_ptr->evflag)
|
||||
if (system->pair_ptr->evflag)
|
||||
system->pair_ptr->ev_tally(i,i,system->n,1,0.0,en_tmp,0.0,0.0,0.0,0.0);
|
||||
}
|
||||
}
|
||||
|
|
|
@ -36,11 +36,11 @@ void Reset_Atoms( reax_system* system, control_params *control )
|
|||
reax_atom *atom;
|
||||
|
||||
system->numH = 0;
|
||||
if( control->hbond_cut > 0 )
|
||||
if (control->hbond_cut > 0)
|
||||
for( i = 0; i < system->n; ++i ) {
|
||||
atom = &(system->my_atoms[i]);
|
||||
if (atom->type < 0) continue;
|
||||
if( system->reax_param.sbp[ atom->type ].p_hbond == 1 )
|
||||
if (system->reax_param.sbp[ atom->type ].p_hbond == 1)
|
||||
atom->Hindex = system->numH++;
|
||||
else atom->Hindex = -1;
|
||||
}
|
||||
|
@ -139,9 +139,9 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* is reallocation needed? */
|
||||
if( total_bonds >= bonds->num_intrs * DANGER_ZONE ) {
|
||||
if (total_bonds >= bonds->num_intrs * DANGER_ZONE) {
|
||||
workspace->realloc.bonds = 1;
|
||||
if( total_bonds >= bonds->num_intrs ) {
|
||||
if (total_bonds >= bonds->num_intrs) {
|
||||
fprintf(stderr,
|
||||
"p%d: not enough space for bonds! total=%d allocated=%d\n",
|
||||
system->my_rank, total_bonds, bonds->num_intrs );
|
||||
|
@ -150,14 +150,14 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
|
|||
}
|
||||
}
|
||||
|
||||
if( control->hbond_cut > 0 && system->numH > 0 ) {
|
||||
if (control->hbond_cut > 0 && system->numH > 0) {
|
||||
hbonds = (*lists) + HBONDS;
|
||||
total_hbonds = 0;
|
||||
|
||||
/* reset start-end indexes */
|
||||
for( i = 0; i < system->n; ++i ) {
|
||||
Hindex = system->my_atoms[i].Hindex;
|
||||
if( Hindex > -1 ) {
|
||||
if (Hindex > -1) {
|
||||
Set_Start_Index( Hindex, total_hbonds, hbonds );
|
||||
Set_End_Index( Hindex, total_hbonds, hbonds );
|
||||
total_hbonds += system->my_atoms[i].num_hbonds;
|
||||
|
@ -165,9 +165,9 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* is reallocation needed? */
|
||||
if( total_hbonds >= hbonds->num_intrs * 0.90/*DANGER_ZONE*/ ) {
|
||||
if (total_hbonds >= hbonds->num_intrs * 0.90/*DANGER_ZONE*/) {
|
||||
workspace->realloc.hbonds = 1;
|
||||
if( total_hbonds >= hbonds->num_intrs ) {
|
||||
if (total_hbonds >= hbonds->num_intrs) {
|
||||
fprintf(stderr,
|
||||
"p%d: not enough space for hbonds! total=%d allocated=%d\n",
|
||||
system->my_rank, total_hbonds, hbonds->num_intrs );
|
||||
|
|
|
@ -60,7 +60,7 @@ void Compute_System_Energy( reax_system *system, simulation_data *data,
|
|||
|
||||
data->my_en.e_tot = data->my_en.e_pot + E_CONV * data->my_en.e_kin;
|
||||
|
||||
if( system->my_rank == MASTER_NODE ) {
|
||||
if (system->my_rank == MASTER_NODE) {
|
||||
data->sys_en.e_bond = sys_en[0];
|
||||
data->sys_en.e_ov = sys_en[1];
|
||||
data->sys_en.e_un = sys_en[2];
|
||||
|
|
|
@ -58,7 +58,7 @@ void *smalloc( rc_bigint n, const char *name, MPI_Comm comm )
|
|||
{
|
||||
void *ptr;
|
||||
|
||||
if( n <= 0 ) {
|
||||
if (n <= 0) {
|
||||
fprintf( stderr, "WARNING: trying to allocate %ld bytes for array %s. ",
|
||||
n, name );
|
||||
fprintf( stderr, "returning NULL.\n" );
|
||||
|
@ -66,7 +66,7 @@ void *smalloc( rc_bigint n, const char *name, MPI_Comm comm )
|
|||
}
|
||||
|
||||
ptr = malloc( n );
|
||||
if( ptr == NULL ) {
|
||||
if (ptr == NULL) {
|
||||
fprintf( stderr, "ERROR: failed to allocate %ld bytes for array %s",
|
||||
n, name );
|
||||
MPI_Abort( comm, INSUFFICIENT_MEMORY );
|
||||
|
@ -81,14 +81,14 @@ void *scalloc( rc_bigint n, rc_bigint size, const char *name, MPI_Comm comm )
|
|||
{
|
||||
void *ptr;
|
||||
|
||||
if( n <= 0 ) {
|
||||
if (n <= 0) {
|
||||
fprintf( stderr, "WARNING: trying to allocate %ld elements for array %s. ",
|
||||
n, name );
|
||||
fprintf( stderr, "returning NULL.\n" );
|
||||
return NULL;
|
||||
}
|
||||
|
||||
if( size <= 0 ) {
|
||||
if (size <= 0) {
|
||||
fprintf( stderr, "WARNING: elements size for array %s is %ld. ",
|
||||
name, size );
|
||||
fprintf( stderr, "returning NULL.\n" );
|
||||
|
@ -96,7 +96,7 @@ void *scalloc( rc_bigint n, rc_bigint size, const char *name, MPI_Comm comm )
|
|||
}
|
||||
|
||||
ptr = calloc( n, size );
|
||||
if( ptr == NULL ) {
|
||||
if (ptr == NULL) {
|
||||
fprintf( stderr, "ERROR: failed to allocate %ld bytes for array %s",
|
||||
n*size, name );
|
||||
MPI_Abort( comm, INSUFFICIENT_MEMORY );
|
||||
|
@ -109,7 +109,7 @@ void *scalloc( rc_bigint n, rc_bigint size, const char *name, MPI_Comm comm )
|
|||
/* safe free */
|
||||
void sfree( void *ptr, const char *name )
|
||||
{
|
||||
if( ptr == NULL ) {
|
||||
if (ptr == NULL) {
|
||||
fprintf( stderr, "WARNING: trying to free the already NULL pointer %s!\n",
|
||||
name );
|
||||
return;
|
||||
|
|
|
@ -74,19 +74,19 @@ double Calculate_Omega( rvec dvec_ij, double r_ij,
|
|||
hnhe = r_ij * r_kl * sin_ijk * cos_jkl;
|
||||
|
||||
poem = 2.0 * r_ij * r_kl * sin_ijk * sin_jkl;
|
||||
if( poem < 1e-20 ) poem = 1e-20;
|
||||
if (poem < 1e-20) poem = 1e-20;
|
||||
|
||||
tel = SQR( r_ij ) + SQR( r_jk ) + SQR( r_kl ) - SQR( r_li ) -
|
||||
2.0 * ( r_ij * r_jk * cos_ijk - r_ij * r_kl * cos_ijk * cos_jkl +
|
||||
r_jk * r_kl * cos_jkl );
|
||||
|
||||
arg = tel / poem;
|
||||
if( arg > 1.0 ) arg = 1.0;
|
||||
if( arg < -1.0 ) arg = -1.0;
|
||||
if (arg > 1.0) arg = 1.0;
|
||||
if (arg < -1.0) arg = -1.0;
|
||||
|
||||
if( sin_ijk >= 0 && sin_ijk <= MIN_SINE ) sin_ijk = MIN_SINE;
|
||||
if (sin_ijk >= 0 && sin_ijk <= MIN_SINE) sin_ijk = MIN_SINE;
|
||||
else if( sin_ijk <= 0 && sin_ijk >= -MIN_SINE ) sin_ijk = -MIN_SINE;
|
||||
if( sin_jkl >= 0 && sin_jkl <= MIN_SINE ) sin_jkl = MIN_SINE;
|
||||
if (sin_jkl >= 0 && sin_jkl <= MIN_SINE) sin_jkl = MIN_SINE;
|
||||
else if( sin_jkl <= 0 && sin_jkl >= -MIN_SINE ) sin_jkl = -MIN_SINE;
|
||||
|
||||
// dcos_omega_di
|
||||
|
@ -181,9 +181,9 @@ void Torsion_Angles( reax_system *system, control_params *control,
|
|||
bo_jk = &( pbond_jk->bo_data );
|
||||
BOA_jk = bo_jk->BO - control->thb_cut;
|
||||
|
||||
if( system->my_atoms[j].orig_id > system->my_atoms[k].orig_id )
|
||||
if (system->my_atoms[j].orig_id > system->my_atoms[k].orig_id)
|
||||
continue;
|
||||
if( system->my_atoms[j].orig_id == system->my_atoms[k].orig_id ) {
|
||||
if (system->my_atoms[j].orig_id == system->my_atoms[k].orig_id) {
|
||||
if (system->my_atoms[k].x[2] < system->my_atoms[j].x[2]) continue;
|
||||
if (system->my_atoms[k].x[2] == system->my_atoms[j].x[2] &&
|
||||
system->my_atoms[k].x[1] < system->my_atoms[j].x[1]) continue;
|
||||
|
@ -192,10 +192,10 @@ void Torsion_Angles( reax_system *system, control_params *control,
|
|||
system->my_atoms[k].x[0] < system->my_atoms[j].x[0]) continue;
|
||||
}
|
||||
|
||||
if( bo_jk->BO > control->thb_cut/*0*/ && Num_Entries(pk, thb_intrs) ) {
|
||||
if (bo_jk->BO > control->thb_cut/*0*/ && Num_Entries(pk, thb_intrs)) {
|
||||
pj = pbond_jk->sym_index; // pj points to j on k's list
|
||||
|
||||
if( Num_Entries(pj, thb_intrs) ) {
|
||||
if (Num_Entries(pj, thb_intrs)) {
|
||||
type_k = system->my_atoms[k].type;
|
||||
Delta_k = workspace->Delta_boc[k];
|
||||
r_jk = pbond_jk->d;
|
||||
|
@ -218,7 +218,7 @@ void Torsion_Angles( reax_system *system, control_params *control,
|
|||
pbond_ij = &( bonds->select.bond_list[pij] );
|
||||
bo_ij = &( pbond_ij->bo_data );
|
||||
|
||||
if( bo_ij->BO > control->thb_cut/*0*/ ) {
|
||||
if (bo_ij->BO > control->thb_cut/*0*/) {
|
||||
i = p_ijk->thb;
|
||||
type_i = system->my_atoms[i].type;
|
||||
r_ij = pbond_ij->d;
|
||||
|
@ -228,7 +228,7 @@ void Torsion_Angles( reax_system *system, control_params *control,
|
|||
sin_ijk = sin( theta_ijk );
|
||||
cos_ijk = cos( theta_ijk );
|
||||
//tan_ijk_i = 1. / tan( theta_ijk );
|
||||
if( sin_ijk >= 0 && sin_ijk <= MIN_SINE )
|
||||
if (sin_ijk >= 0 && sin_ijk <= MIN_SINE)
|
||||
tan_ijk_i = cos_ijk / MIN_SINE;
|
||||
else if( sin_ijk <= 0 && sin_ijk >= -MIN_SINE )
|
||||
tan_ijk_i = cos_ijk / -MIN_SINE;
|
||||
|
@ -260,7 +260,7 @@ void Torsion_Angles( reax_system *system, control_params *control,
|
|||
sin_jkl = sin( theta_jkl );
|
||||
cos_jkl = cos( theta_jkl );
|
||||
//tan_jkl_i = 1. / tan( theta_jkl );
|
||||
if( sin_jkl >= 0 && sin_jkl <= MIN_SINE )
|
||||
if (sin_jkl >= 0 && sin_jkl <= MIN_SINE)
|
||||
tan_jkl_i = cos_jkl / MIN_SINE;
|
||||
else if( sin_jkl <= 0 && sin_jkl >= -MIN_SINE )
|
||||
tan_jkl_i = cos_jkl / -MIN_SINE;
|
||||
|
@ -357,7 +357,7 @@ void Torsion_Angles( reax_system *system, control_params *control,
|
|||
bo_jk->Cdbo += (CEtors5 + CEconj2);
|
||||
bo_kl->Cdbo += (CEtors6 + CEconj3);
|
||||
|
||||
if( control->virial == 0 ) {
|
||||
if (control->virial == 0) {
|
||||
/* dcos_theta_ijk */
|
||||
rvec_ScaledAdd( workspace->f[i],
|
||||
CEtors7 + CEconj4, p_ijk->dcos_dk );
|
||||
|
@ -438,7 +438,7 @@ void Torsion_Angles( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* tally into per-atom virials */
|
||||
if( system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
|
||||
if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
|
||||
|
||||
// acquire vectors
|
||||
rvec_ScaledSum( delil, 1., system->my_atoms[l].x,
|
||||
|
@ -463,9 +463,9 @@ void Torsion_Angles( reax_system *system, control_params *control,
|
|||
|
||||
// tally
|
||||
eng_tmp = e_tor + e_con;
|
||||
if( system->pair_ptr->evflag)
|
||||
if (system->pair_ptr->evflag)
|
||||
system->pair_ptr->ev_tally(j,k,natoms,1,eng_tmp,0.0,0.0,0.0,0.0,0.0);
|
||||
if( system->pair_ptr->vflag_atom)
|
||||
if (system->pair_ptr->vflag_atom)
|
||||
system->pair_ptr->v_tally4(i,j,k,l,fi_tmp,fj_tmp,fk_tmp,delil,deljl,delkl);
|
||||
}
|
||||
} // pl check ends
|
||||
|
|
|
@ -32,12 +32,12 @@
|
|||
int Reallocate_Output_Buffer( output_controls *out_control, int req_space,
|
||||
MPI_Comm comm )
|
||||
{
|
||||
if( out_control->buffer_len > 0 )
|
||||
if (out_control->buffer_len > 0)
|
||||
free( out_control->buffer );
|
||||
|
||||
out_control->buffer_len = (int)(req_space*SAFE_ZONE);
|
||||
out_control->buffer = (char*) malloc(out_control->buffer_len*sizeof(char));
|
||||
if( out_control->buffer == NULL ) {
|
||||
if (out_control->buffer == NULL) {
|
||||
fprintf( stderr,
|
||||
"insufficient memory for required buffer size %d. terminating!\n",
|
||||
(int) (req_space*SAFE_ZONE) );
|
||||
|
@ -51,7 +51,7 @@ int Reallocate_Output_Buffer( output_controls *out_control, int req_space,
|
|||
void Write_Skip_Line( output_controls *out_control, mpi_datatypes * /*mpi_data*/,
|
||||
int my_rank, int skip, int num_section )
|
||||
{
|
||||
if( my_rank == MASTER_NODE )
|
||||
if (my_rank == MASTER_NODE)
|
||||
fprintf( out_control->strj, INT2_LINE,
|
||||
"chars_to_skip_section:", skip, num_section );
|
||||
}
|
||||
|
@ -82,11 +82,11 @@ int Write_Header( reax_system *system, control_params *control,
|
|||
num_hdr_lines = NUM_HEADER_LINES;
|
||||
my_hdr_lines = num_hdr_lines * ( system->my_rank == MASTER_NODE );
|
||||
buffer_req = my_hdr_lines * HEADER_LINE_LEN;
|
||||
if( buffer_req > out_control->buffer_len * DANGER_ZONE )
|
||||
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
|
||||
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
|
||||
|
||||
/* only the master node writes into trajectory header */
|
||||
if( system->my_rank == MASTER_NODE ) {
|
||||
if (system->my_rank == MASTER_NODE) {
|
||||
/* clear the contents of line & buffer */
|
||||
out_control->line[0] = 0;
|
||||
out_control->buffer[0] = 0;
|
||||
|
@ -252,7 +252,7 @@ int Write_Header( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* dump out the buffer */
|
||||
if( system->my_rank == MASTER_NODE )
|
||||
if (system->my_rank == MASTER_NODE)
|
||||
fprintf( out_control->strj, "%s", out_control->buffer );
|
||||
|
||||
return SUCCESS;
|
||||
|
@ -273,11 +273,11 @@ int Write_Init_Desc( reax_system *system, control_params * /*control*/,
|
|||
Write_Skip_Line( out_control, mpi_data, me,
|
||||
system->bigN * INIT_DESC_LEN, system->bigN );
|
||||
|
||||
if( out_control->traj_method == REG_TRAJ && me == MASTER_NODE )
|
||||
if (out_control->traj_method == REG_TRAJ && me == MASTER_NODE)
|
||||
buffer_req = system->bigN * INIT_DESC_LEN + 1;
|
||||
else buffer_req = system->n * INIT_DESC_LEN + 1;
|
||||
|
||||
if( buffer_req > out_control->buffer_len * DANGER_ZONE )
|
||||
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
|
||||
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
|
||||
|
||||
out_control->line[0] = 0;
|
||||
|
@ -297,7 +297,7 @@ int Write_Init_Desc( reax_system *system, control_params * /*control*/,
|
|||
} else {
|
||||
buffer_len = system->n * INIT_DESC_LEN;
|
||||
for( i = 0; i < np; ++i )
|
||||
if( i != MASTER_NODE ) {
|
||||
if (i != MASTER_NODE) {
|
||||
MPI_Recv( out_control->buffer + buffer_len, buffer_req - buffer_len,
|
||||
MPI_CHAR, i, np*INIT_DESCS+i, mpi_data->world, &status );
|
||||
MPI_Get_count( &status, MPI_CHAR, &cnt );
|
||||
|
@ -341,8 +341,8 @@ int Init_Traj( reax_system *system, control_params *control,
|
|||
out_control->buffer = NULL;
|
||||
|
||||
/* write trajectory header and atom info, if applicable */
|
||||
if( out_control->traj_method == REG_TRAJ) {
|
||||
if( system->my_rank == MASTER_NODE )
|
||||
if (out_control->traj_method == REG_TRAJ) {
|
||||
if (system->my_rank == MASTER_NODE)
|
||||
out_control->strj = fopen( fname, "w" );
|
||||
} else {
|
||||
strcpy( msg, "init_traj: unknown trajectory option" );
|
||||
|
@ -366,11 +366,11 @@ int Write_Frame_Header( reax_system *system, control_params *control,
|
|||
num_frm_hdr_lines = 22;
|
||||
my_frm_hdr_lines = num_frm_hdr_lines * ( me == MASTER_NODE );
|
||||
buffer_req = my_frm_hdr_lines * HEADER_LINE_LEN;
|
||||
if( buffer_req > out_control->buffer_len * DANGER_ZONE )
|
||||
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
|
||||
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
|
||||
|
||||
/* only the master node writes into trajectory header */
|
||||
if( me == MASTER_NODE ) {
|
||||
if (me == MASTER_NODE) {
|
||||
/* clear the contents of line & buffer */
|
||||
out_control->line[0] = 0;
|
||||
out_control->buffer[0] = 0;
|
||||
|
@ -473,7 +473,7 @@ int Write_Frame_Header( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* dump out the buffer */
|
||||
if( system->my_rank == MASTER_NODE )
|
||||
if (system->my_rank == MASTER_NODE)
|
||||
fprintf( out_control->strj, "%s", out_control->buffer );
|
||||
|
||||
return SUCCESS;
|
||||
|
@ -495,11 +495,11 @@ int Write_Atoms( reax_system *system, control_params * /*control*/,
|
|||
Write_Skip_Line( out_control, mpi_data, me,
|
||||
system->bigN*line_len, system->bigN );
|
||||
|
||||
if( out_control->traj_method == REG_TRAJ && me == MASTER_NODE )
|
||||
if (out_control->traj_method == REG_TRAJ && me == MASTER_NODE)
|
||||
buffer_req = system->bigN * line_len + 1;
|
||||
else buffer_req = system->n * line_len + 1;
|
||||
|
||||
if( buffer_req > out_control->buffer_len * DANGER_ZONE )
|
||||
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
|
||||
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
|
||||
|
||||
/* fill in buffer */
|
||||
|
@ -545,7 +545,7 @@ int Write_Atoms( reax_system *system, control_params * /*control*/,
|
|||
} else {
|
||||
buffer_len = system->n * line_len;
|
||||
for( i = 0; i < np; ++i )
|
||||
if( i != MASTER_NODE ) {
|
||||
if (i != MASTER_NODE) {
|
||||
MPI_Recv( out_control->buffer + buffer_len, buffer_req - buffer_len,
|
||||
MPI_CHAR, i, np*ATOM_LINES+i, mpi_data->world, &status );
|
||||
MPI_Get_count( &status, MPI_CHAR, &cnt );
|
||||
|
@ -587,11 +587,11 @@ int Write_Bonds(reax_system *system, control_params *control, reax_list *bonds,
|
|||
|
||||
Write_Skip_Line( out_control, mpi_data, me, num_bonds*line_len, num_bonds );
|
||||
|
||||
if( out_control->traj_method == REG_TRAJ && me == MASTER_NODE )
|
||||
if (out_control->traj_method == REG_TRAJ && me == MASTER_NODE)
|
||||
buffer_req = num_bonds * line_len + 1;
|
||||
else buffer_req = my_bonds * line_len + 1;
|
||||
|
||||
if( buffer_req > out_control->buffer_len * DANGER_ZONE )
|
||||
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
|
||||
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
|
||||
|
||||
/* fill in the buffer */
|
||||
|
@ -635,7 +635,7 @@ int Write_Bonds(reax_system *system, control_params *control, reax_list *bonds,
|
|||
} else {
|
||||
buffer_len = my_bonds * line_len;
|
||||
for( i = 0; i < np; ++i )
|
||||
if( i != MASTER_NODE ) {
|
||||
if (i != MASTER_NODE) {
|
||||
MPI_Recv( out_control->buffer + buffer_len, buffer_req - buffer_len,
|
||||
MPI_CHAR, i, np*BOND_LINES+i, mpi_data->world, &status );
|
||||
MPI_Get_count( &status, MPI_CHAR, &cnt );
|
||||
|
@ -671,7 +671,7 @@ int Write_Angles( reax_system *system, control_params *control,
|
|||
bo_ij = &(bonds->select.bond_list[pi]);
|
||||
i = bo_ij->nbr;
|
||||
|
||||
if( bo_ij->bo_data.BO >= control->bg_cut ) // physical j&i bond
|
||||
if (bo_ij->bo_data.BO >= control->bg_cut) // physical j&i bond
|
||||
for( pk = Start_Index( pi, thb_intrs );
|
||||
pk < End_Index( pi, thb_intrs ); ++pk ) {
|
||||
angle_ijk = &(thb_intrs->select.three_body_list[pk]);
|
||||
|
@ -688,11 +688,11 @@ int Write_Angles( reax_system *system, control_params *control,
|
|||
|
||||
Write_Skip_Line( out_control, mpi_data, me, num_angles*line_len, num_angles );
|
||||
|
||||
if( out_control->traj_method == REG_TRAJ && me == MASTER_NODE )
|
||||
if (out_control->traj_method == REG_TRAJ && me == MASTER_NODE)
|
||||
buffer_req = num_angles * line_len + 1;
|
||||
else buffer_req = my_angles * line_len + 1;
|
||||
|
||||
if( buffer_req > out_control->buffer_len * DANGER_ZONE )
|
||||
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
|
||||
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
|
||||
|
||||
/* fill in the buffer */
|
||||
|
@ -704,7 +704,7 @@ int Write_Angles( reax_system *system, control_params *control,
|
|||
bo_ij = &(bonds->select.bond_list[pi]);
|
||||
i = bo_ij->nbr;
|
||||
|
||||
if( bo_ij->bo_data.BO >= control->bg_cut ) // physical j&i bond
|
||||
if (bo_ij->bo_data.BO >= control->bg_cut) // physical j&i bond
|
||||
for( pk = Start_Index( pi, thb_intrs );
|
||||
pk < End_Index( pi, thb_intrs ); ++pk ) {
|
||||
angle_ijk = &(thb_intrs->select.three_body_list[pk]);
|
||||
|
@ -730,7 +730,7 @@ int Write_Angles( reax_system *system, control_params *control,
|
|||
} else {
|
||||
buffer_len = my_angles * line_len;
|
||||
for( i = 0; i < np; ++i )
|
||||
if( i != MASTER_NODE ) {
|
||||
if (i != MASTER_NODE) {
|
||||
MPI_Recv( out_control->buffer + buffer_len, buffer_req - buffer_len,
|
||||
MPI_CHAR, i, np*ANGLE_LINES+i, mpi_data->world, &status );
|
||||
MPI_Get_count( &status, MPI_CHAR, &cnt );
|
||||
|
@ -750,13 +750,13 @@ int Append_Frame( reax_system *system, control_params *control,
|
|||
{
|
||||
Write_Frame_Header( system, control, data, out_control, mpi_data );
|
||||
|
||||
if( out_control->write_atoms )
|
||||
if (out_control->write_atoms)
|
||||
Write_Atoms( system, control, out_control, mpi_data );
|
||||
|
||||
if( out_control->write_bonds )
|
||||
if (out_control->write_bonds)
|
||||
Write_Bonds( system, control, (*lists + BONDS), out_control, mpi_data );
|
||||
|
||||
if( out_control->write_angles )
|
||||
if (out_control->write_angles)
|
||||
Write_Angles( system, control, (*lists + BONDS), (*lists + THREE_BODIES),
|
||||
out_control, mpi_data );
|
||||
|
||||
|
@ -766,7 +766,7 @@ int Append_Frame( reax_system *system, control_params *control,
|
|||
|
||||
int End_Traj( int my_rank, output_controls *out_control )
|
||||
{
|
||||
if( my_rank == MASTER_NODE )
|
||||
if (my_rank == MASTER_NODE)
|
||||
fclose( out_control->strj );
|
||||
|
||||
free( out_control->buffer );
|
||||
|
|
|
@ -44,8 +44,8 @@ void Calculate_Theta( rvec dvec_ji, double d_ji, rvec dvec_jk, double d_jk,
|
|||
double *theta, double *cos_theta )
|
||||
{
|
||||
(*cos_theta) = Dot( dvec_ji, dvec_jk, 3 ) / ( d_ji * d_jk );
|
||||
if( *cos_theta > 1. ) *cos_theta = 1.0;
|
||||
if( *cos_theta < -1. ) *cos_theta = -1.0;
|
||||
if (*cos_theta > 1.) *cos_theta = 1.0;
|
||||
if (*cos_theta < -1.) *cos_theta = -1.0;
|
||||
|
||||
(*theta) = acos( *cos_theta );
|
||||
}
|
||||
|
@ -151,13 +151,13 @@ void Valence_Angles( reax_system *system, control_params *control,
|
|||
SBO = SBOp + (1 - prod_SBO) * (-workspace->Delta_boc[j] - p_val8 * vlpadj);
|
||||
dSBO1 = -8 * prod_SBO * ( workspace->Delta_boc[j] + p_val8 * vlpadj );
|
||||
|
||||
if( SBO <= 0 )
|
||||
if (SBO <= 0)
|
||||
SBO2 = 0, CSBO2 = 0;
|
||||
else if( SBO > 0 && SBO <= 1 ) {
|
||||
else if (SBO > 0 && SBO <= 1) {
|
||||
SBO2 = pow( SBO, p_val9 );
|
||||
CSBO2 = p_val9 * pow( SBO, p_val9 - 1 );
|
||||
}
|
||||
else if( SBO > 1 && SBO < 2 ) {
|
||||
else if (SBO > 1 && SBO < 2) {
|
||||
SBO2 = 2 - pow( 2-SBO, p_val9 );
|
||||
CSBO2 = p_val9 * pow( 2 - SBO, p_val9 - 1 );
|
||||
}
|
||||
|
@ -183,7 +183,7 @@ void Valence_Angles( reax_system *system, control_params *control,
|
|||
end_pk = End_Index( pk, thb_intrs );
|
||||
|
||||
for( t = start_pk; t < end_pk; ++t )
|
||||
if( thb_intrs->select.three_body_list[t].thb == i ) {
|
||||
if (thb_intrs->select.three_body_list[t].thb == i) {
|
||||
p_ijk = &(thb_intrs->select.three_body_list[num_thb_intrs] );
|
||||
p_kji = &(thb_intrs->select.three_body_list[t]);
|
||||
|
||||
|
@ -220,20 +220,20 @@ void Valence_Angles( reax_system *system, control_params *control,
|
|||
p_ijk->theta = theta;
|
||||
|
||||
sin_theta = sin( theta );
|
||||
if( sin_theta < 1.0e-5 )
|
||||
if (sin_theta < 1.0e-5)
|
||||
sin_theta = 1.0e-5;
|
||||
|
||||
++num_thb_intrs;
|
||||
|
||||
|
||||
if( (j < system->n) && (BOA_jk > 0.0) &&
|
||||
if ((j < system->n) && (BOA_jk > 0.0) &&
|
||||
(bo_ij->BO > control->thb_cut) &&
|
||||
(bo_jk->BO > control->thb_cut) &&
|
||||
(bo_ij->BO * bo_jk->BO > control->thb_cutsq) ) {
|
||||
thbh = &( system->reax_param.thbp[ type_i ][ type_j ][ type_k ] );
|
||||
|
||||
for( cnt = 0; cnt < thbh->cnt; ++cnt ) {
|
||||
if( fabs(thbh->prm[cnt].p_val1) > 0.001 ) {
|
||||
if (fabs(thbh->prm[cnt].p_val1) > 0.001) {
|
||||
thbp = &( thbh->prm[cnt] );
|
||||
|
||||
/* ANGLE ENERGY */
|
||||
|
@ -263,7 +263,7 @@ void Valence_Angles( reax_system *system, control_params *control,
|
|||
theta_0 = DEG2RAD( theta_0 );
|
||||
|
||||
expval2theta = exp( -p_val2 * SQR(theta_0 - theta) );
|
||||
if( p_val1 >= 0 )
|
||||
if (p_val1 >= 0)
|
||||
expval12theta = p_val1 * (1.0 - expval2theta);
|
||||
else // To avoid linear Me-H-Me angles (6/6/06)
|
||||
expval12theta = p_val1 * -expval2theta;
|
||||
|
@ -354,7 +354,7 @@ void Valence_Angles( reax_system *system, control_params *control,
|
|||
bo_jt->Cdbopi2 += CEval5;
|
||||
}
|
||||
|
||||
if( control->virial == 0 ) {
|
||||
if (control->virial == 0) {
|
||||
rvec_ScaledAdd( workspace->f[i], CEval8, p_ijk->dcos_di );
|
||||
rvec_ScaledAdd( workspace->f[j], CEval8, p_ijk->dcos_dj );
|
||||
rvec_ScaledAdd( workspace->f[k], CEval8, p_ijk->dcos_dk );
|
||||
|
@ -373,7 +373,7 @@ void Valence_Angles( reax_system *system, control_params *control,
|
|||
}
|
||||
|
||||
/* tally into per-atom virials */
|
||||
if( system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
|
||||
if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
|
||||
|
||||
/* Acquire vectors */
|
||||
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
|
||||
|
@ -387,9 +387,9 @@ void Valence_Angles( reax_system *system, control_params *control,
|
|||
|
||||
eng_tmp = e_ang + e_pen + e_coa;
|
||||
|
||||
if( system->pair_ptr->evflag)
|
||||
if (system->pair_ptr->evflag)
|
||||
system->pair_ptr->ev_tally(j,j,system->N,1,eng_tmp,0.0,0.0,0.0,0.0,0.0);
|
||||
if( system->pair_ptr->vflag_atom)
|
||||
if (system->pair_ptr->vflag_atom)
|
||||
system->pair_ptr->v_tally3(i,j,k,fi_tmp,fk_tmp,delij,delkj);
|
||||
}
|
||||
}
|
||||
|
@ -402,9 +402,9 @@ void Valence_Angles( reax_system *system, control_params *control,
|
|||
}
|
||||
}
|
||||
|
||||
if( num_thb_intrs >= thb_intrs->num_intrs * DANGER_ZONE ) {
|
||||
if (num_thb_intrs >= thb_intrs->num_intrs * DANGER_ZONE) {
|
||||
workspace->realloc.num_3body = num_thb_intrs;
|
||||
if( num_thb_intrs > thb_intrs->num_intrs ) {
|
||||
if (num_thb_intrs > thb_intrs->num_intrs) {
|
||||
fprintf( stderr, "step%d-ran out of space on angle_list: top=%d, max=%d",
|
||||
data->step, num_thb_intrs, thb_intrs->num_intrs );
|
||||
MPI_Abort( MPI_COMM_WORLD, INSUFFICIENT_MEMORY );
|
||||
|
|
|
@ -105,7 +105,7 @@ void rtensor_MatVec( rvec ret, rtensor m, rvec v )
|
|||
int i;
|
||||
rvec temp;
|
||||
|
||||
if( ret == v )
|
||||
if (ret == v)
|
||||
{
|
||||
for( i = 0; i < 3; ++i )
|
||||
temp[i] = m[i][0] * v[0] + m[i][1] * v[1] + m[i][2] * v[2];
|
||||
|
|
|
@ -90,7 +90,7 @@ ComputeVoronoi::ComputeVoronoi(LAMMPS *lmp, int narg, char **arg) :
|
|||
else if (strcmp(arg[iarg], "surface") == 0) {
|
||||
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
|
||||
// group all is a special case where we just skip group testing
|
||||
if(strcmp(arg[iarg+1], "all") == 0) {
|
||||
if (strcmp(arg[iarg+1], "all") == 0) {
|
||||
surface = VOROSURF_ALL;
|
||||
} else {
|
||||
sgroup = group->find(arg[iarg+1]);
|
||||
|
@ -266,7 +266,7 @@ void ComputeVoronoi::buildCells()
|
|||
double **x = atom->x;
|
||||
|
||||
// setup bounds for voro++ domain for orthogonal and triclinic simulation boxes
|
||||
if( domain->triclinic ) {
|
||||
if (domain->triclinic) {
|
||||
// triclinic box: embed parallelepiped into orthogonal voro++ domain
|
||||
|
||||
// cutghost is in lamda coordinates for triclinic boxes, use subxx_lamda
|
||||
|
@ -356,7 +356,7 @@ void ComputeVoronoi::buildCells()
|
|||
|
||||
// pass coordinates for local and ghost atoms to voro++
|
||||
for (i = 0; i < nall; i++) {
|
||||
if( !onlyGroup || (mask[i] & groupbit) )
|
||||
if (!onlyGroup || (mask[i] & groupbit))
|
||||
con_poly->put(i,x[i][0],x[i][1],x[i][2],rfield[i]);
|
||||
}
|
||||
} else {
|
||||
|
@ -374,7 +374,7 @@ void ComputeVoronoi::buildCells()
|
|||
|
||||
// pass coordinates for local and ghost atoms to voro++
|
||||
for (i = 0; i < nall; i++)
|
||||
if( !onlyGroup || (mask[i] & groupbit) )
|
||||
if (!onlyGroup || (mask[i] & groupbit))
|
||||
con_mono->put(i,x[i][0],x[i][1],x[i][2]);
|
||||
}
|
||||
}
|
||||
|
@ -625,7 +625,7 @@ double ComputeVoronoi::memory_usage()
|
|||
void ComputeVoronoi::compute_vector()
|
||||
{
|
||||
invoked_vector = update->ntimestep;
|
||||
if( invoked_peratom < invoked_vector ) compute_peratom();
|
||||
if (invoked_peratom < invoked_vector) compute_peratom();
|
||||
|
||||
for( int i=0; i<size_vector; ++i ) sendvector[i] = edge[i];
|
||||
MPI_Allreduce(sendvector,edge,size_vector,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
@ -636,7 +636,7 @@ void ComputeVoronoi::compute_vector()
|
|||
void ComputeVoronoi::compute_local()
|
||||
{
|
||||
invoked_local = update->ntimestep;
|
||||
if( invoked_peratom < invoked_local ) compute_peratom();
|
||||
if (invoked_peratom < invoked_local) compute_peratom();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
Loading…
Reference in New Issue