Skip to content
Snippets Groups Projects
Commit 20ece267 authored by Kurt A. O'Hearn's avatar Kurt A. O'Hearn
Browse files

PG-PuReMD: backport missing valence angle correction from sPuReMD.

parent 5f6a651b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
PG-PuReMD/src/cuda/cuda_valence_angles.cu
+ 1
− 1
View file @ 20ece267
......@@ -114,7 +114,7 @@ CUDA_GLOBAL void Cuda_Valence_Angles_Part1( reax_atom *my_atoms,
else
{
vlpadj = workspace.nlp[j];
dSBO2 = (prod_SBO - 1.0) * (1.0 - p_val8 * workspace.dDelta_lp[j]);
dSBO2 = (prod_SBO - 1.0) * (1.0 + p_val8 * workspace.dDelta_lp[j]);
}
SBO = SBOp + (1.0 - prod_SBO) * (-workspace.Delta_boc[j] - p_val8 * vlpadj);
......
PG-PuReMD/src/nonbonded.c
+ 1
− 1
View file @ 20ece267
......@@ -206,7 +206,7 @@ void vdW_Coulomb_Energy( reax_system const * const system,
if ( control->virial == 0 )
{
rvec_ScaledAdd( workspace->f[i],
-(CEvd + CEclmb) / r_ij, far_nbr_list->far_nbr_list.dvec[pj] );
-1.0 * (CEvd + CEclmb) / r_ij, far_nbr_list->far_nbr_list.dvec[pj] );
rvec_ScaledAdd( workspace->f[j],
(CEvd + CEclmb) / r_ij, far_nbr_list->far_nbr_list.dvec[pj] );
}
......
PG-PuReMD/src/valence_angles.c
+ 244
− 242
View file @ 20ece267
......@@ -186,7 +186,7 @@ void Valence_Angles( reax_system * const system, control_params * const control,
else
{
vlpadj = workspace->nlp[j];
dSBO2 = (prod_SBO - 1.0) * (1.0 - p_val8 * workspace->dDelta_lp[j]);
dSBO2 = (prod_SBO - 1.0) * (1.0 + p_val8 * workspace->dDelta_lp[j]);
}
SBO = SBOp + (1.0 - prod_SBO) * (-workspace->Delta_boc[j] - p_val8 * vlpadj);
......@@ -300,259 +300,261 @@ void Valence_Angles( reax_system * const system, control_params * const control,
/* Fortran ReaxFF code hard-codes the constant below
* as of 2019-02-27, so use that for now */
if ( j < system->n && BOA_jk >= 0.0 && (bo_ij->BO * bo_jk->BO) >= 0.00001 )
// if ( j < system->n && BOA_jk >= 0.0 && (bo_ij->BO * bo_jk->BO) > SQR(control->thb_cut) )
if ( j >= system->n || BOA_jk < 0.0 || (bo_ij->BO * bo_jk->BO) < 0.00001 )
// if ( j < system->n || BOA_jk < 0.0 || (bo_ij->BO * bo_jk->BO) < SQR(control->thb_cut) )
{
thbh = &system->reax_param.thbp[
index_thbp( type_i, type_j, type_k, system->reax_param.num_atom_types ) ];
continue;
}
thbh = &system->reax_param.thbp[
index_thbp( type_i, type_j, type_k, system->reax_param.num_atom_types ) ];
for ( cnt = 0; cnt < thbh->cnt; ++cnt )
{
/* valence angle does not exist in the force field */
if ( FABS(thbh->prm[cnt].p_val1) < 0.001 )
{
continue;
}
thbp = &thbh->prm[cnt];
/* calculate valence angle energy */
p_val1 = thbp->p_val1;
p_val2 = thbp->p_val2;
p_val4 = thbp->p_val4;
p_val7 = thbp->p_val7;
theta_00 = thbp->theta_00;
exp3ij = EXP( -p_val3 * POW( BOA_ij, p_val4 ) );
f7_ij = 1.0 - exp3ij;
Cf7ij = p_val3 * p_val4
* POW( BOA_ij, p_val4 - 1.0 ) * exp3ij;
exp3jk = EXP( -p_val3 * POW( BOA_jk, p_val4 ) );
f7_jk = 1.0 - exp3jk;
Cf7jk = p_val3 * p_val4
* POW( BOA_jk, p_val4 - 1.0 ) * exp3jk;
expval7 = EXP( -p_val7 * workspace->Delta_boc[j] );
trm8 = 1.0 + expval6 + expval7;
f8_Dj = p_val5 - (p_val5 - 1.0) * (2.0 + expval6) / trm8;
Cf8j = ( (1.0 - p_val5) / SQR(trm8) )
* (p_val6 * expval6 * trm8
- (2.0 + expval6) * ( p_val6 * expval6 - p_val7 * expval7 ));
theta_0 = 180.0 - theta_00 * (1.0 - EXP(-p_val10 * (2.0 - SBO2)));
theta_0 = DEG2RAD( theta_0 );
expval2theta = p_val1 * EXP(-p_val2 * SQR(theta_0 - theta));
if ( p_val1 >= 0.0 )
{
expval12theta = p_val1 - expval2theta;
}
/* To avoid linear Me-H-Me angles (6/6/06) */
else
{
expval12theta = -expval2theta;
}
for ( cnt = 0; cnt < thbh->cnt; ++cnt )
CEval1 = Cf7ij * f7_jk * f8_Dj * expval12theta;
CEval2 = Cf7jk * f7_ij * f8_Dj * expval12theta;
CEval3 = Cf8j * f7_ij * f7_jk * expval12theta;
CEval4 = 2.0 * p_val2 * f7_ij * f7_jk * f8_Dj
* expval2theta * (theta_0 - theta);
Ctheta_0 = p_val10 * DEG2RAD(theta_00)
* EXP( -p_val10 * (2.0 - SBO2) );
CEval5 = CEval4 * Ctheta_0 * CSBO2;
CEval6 = CEval5 * dSBO1;
CEval7 = CEval5 * dSBO2;
CEval8 = CEval4 / sin_theta;
e_ang = f7_ij * f7_jk * f8_Dj * expval12theta;
data->my_en.e_ang += e_ang;
/* calculate penalty for double bonds in valency angles */
p_pen1 = thbp->p_pen1;
exp_pen2ij = EXP( -p_pen2 * SQR( BOA_ij - 2.0 ) );
exp_pen2jk = EXP( -p_pen2 * SQR( BOA_jk - 2.0 ) );
exp_pen3 = EXP( -p_pen3 * workspace->Delta[j] );
exp_pen4 = EXP( p_pen4 * workspace->Delta[j] );
trm_pen34 = 1.0 + exp_pen3 + exp_pen4;
f9_Dj = ( 2.0 + exp_pen3 ) / trm_pen34;
Cf9j = (-p_pen3 * exp_pen3 * trm_pen34
- (2.0 + exp_pen3) * ( -p_pen3 * exp_pen3
+ p_pen4 * exp_pen4 )) / SQR( trm_pen34 );
e_pen = p_pen1 * f9_Dj * exp_pen2ij * exp_pen2jk;
data->my_en.e_pen += e_pen;
CEpen1 = e_pen * Cf9j / f9_Dj;
temp = -2.0 * p_pen2 * e_pen;
CEpen2 = temp * (BOA_ij - 2.0);
CEpen3 = temp * (BOA_jk - 2.0);
/* calculate valency angle conjugation energy */
p_coa1 = thbp->p_coa1;
exp_coa2 = EXP( p_coa2 * workspace->Delta_boc[j] );
e_coa = p_coa1
* EXP( -p_coa4 * SQR(BOA_ij - 1.5) )
* EXP( -p_coa4 * SQR(BOA_jk - 1.5) )
* EXP( -p_coa3 * SQR(workspace->total_bond_order[i] - BOA_ij) )
* EXP( -p_coa3 * SQR(workspace->total_bond_order[k] - BOA_jk) )
/ (1.0 + exp_coa2);
data->my_en.e_coa += e_coa;
CEcoa1 = -2.0 * p_coa4 * (BOA_ij - 1.5) * e_coa;
CEcoa2 = -2.0 * p_coa4 * (BOA_jk - 1.5) * e_coa;
CEcoa3 = -p_coa2 * exp_coa2 * e_coa / (1.0 + exp_coa2);
CEcoa4 = -2.0 * p_coa3 * (workspace->total_bond_order[i] - BOA_ij) * e_coa;
CEcoa5 = -2.0 * p_coa3 * (workspace->total_bond_order[k] - BOA_jk) * e_coa;
/* calculate force contributions */
bo_ij->Cdbo += CEval1 + CEpen2 + (CEcoa1 - CEcoa4);
bo_jk->Cdbo += CEval2 + CEpen3 + (CEcoa2 - CEcoa5);
workspace->CdDelta[j] += (CEval3 + CEval7) + CEpen1 + CEcoa3;
workspace->CdDelta[i] += CEcoa4;
workspace->CdDelta[k] += CEcoa5;
for ( t = start_j; t < end_j; ++t )
{
/* valence angle does not exist in the force field */
if ( FABS(thbh->prm[cnt].p_val1) < 0.001 )
{
continue;
}
thbp = &thbh->prm[cnt];
/* calculate valence angle energy */
p_val1 = thbp->p_val1;
p_val2 = thbp->p_val2;
p_val4 = thbp->p_val4;
p_val7 = thbp->p_val7;
theta_00 = thbp->theta_00;
exp3ij = EXP( -p_val3 * POW( BOA_ij, p_val4 ) );
f7_ij = 1.0 - exp3ij;
Cf7ij = p_val3 * p_val4
* POW( BOA_ij, p_val4 - 1.0 ) * exp3ij;
exp3jk = EXP( -p_val3 * POW( BOA_jk, p_val4 ) );
f7_jk = 1.0 - exp3jk;
Cf7jk = p_val3 * p_val4
* POW( BOA_jk, p_val4 - 1.0 ) * exp3jk;
expval7 = EXP( -p_val7 * workspace->Delta_boc[j] );
trm8 = 1.0 + expval6 + expval7;
f8_Dj = p_val5 - (p_val5 - 1.0) * (2.0 + expval6) / trm8;
Cf8j = ( (1.0 - p_val5) / SQR(trm8) )
* (p_val6 * expval6 * trm8
- (2.0 + expval6) * ( p_val6 * expval6 - p_val7 * expval7 ));
theta_0 = 180.0 - theta_00 * (1.0 - EXP(-p_val10 * (2.0 - SBO2)));
theta_0 = DEG2RAD( theta_0 );
expval2theta = p_val1 * EXP(-p_val2 * SQR(theta_0 - theta));
if ( p_val1 >= 0.0 )
{
expval12theta = p_val1 - expval2theta;
}
/* To avoid linear Me-H-Me angles (6/6/06) */
else
{
expval12theta = -expval2theta;
}
CEval1 = Cf7ij * f7_jk * f8_Dj * expval12theta;
CEval2 = Cf7jk * f7_ij * f8_Dj * expval12theta;
CEval3 = Cf8j * f7_ij * f7_jk * expval12theta;
CEval4 = 2.0 * p_val2 * f7_ij * f7_jk * f8_Dj
* expval2theta * (theta_0 - theta);
Ctheta_0 = p_val10 * DEG2RAD(theta_00)
* EXP( -p_val10 * (2.0 - SBO2) );
CEval5 = CEval4 * Ctheta_0 * CSBO2;
CEval6 = CEval5 * dSBO1;
CEval7 = CEval5 * dSBO2;
CEval8 = CEval4 / sin_theta;
e_ang = f7_ij * f7_jk * f8_Dj * expval12theta;
data->my_en.e_ang += e_ang;
/* calculate penalty for double bonds in valency angles */
p_pen1 = thbp->p_pen1;
exp_pen2ij = EXP( -p_pen2 * SQR( BOA_ij - 2.0 ) );
exp_pen2jk = EXP( -p_pen2 * SQR( BOA_jk - 2.0 ) );
exp_pen3 = EXP( -p_pen3 * workspace->Delta[j] );
exp_pen4 = EXP( p_pen4 * workspace->Delta[j] );
trm_pen34 = 1.0 + exp_pen3 + exp_pen4;
f9_Dj = ( 2.0 + exp_pen3 ) / trm_pen34;
Cf9j = (-p_pen3 * exp_pen3 * trm_pen34
- (2.0 + exp_pen3) * ( -p_pen3 * exp_pen3
+ p_pen4 * exp_pen4 )) / SQR( trm_pen34 );
e_pen = p_pen1 * f9_Dj * exp_pen2ij * exp_pen2jk;
data->my_en.e_pen += e_pen;
CEpen1 = e_pen * Cf9j / f9_Dj;
temp = -2.0 * p_pen2 * e_pen;
CEpen2 = temp * (BOA_ij - 2.0);
CEpen3 = temp * (BOA_jk - 2.0);
/* calculate valency angle conjugation energy */
p_coa1 = thbp->p_coa1;
exp_coa2 = EXP( p_coa2 * workspace->Delta_boc[j] );
e_coa = p_coa1
* EXP( -p_coa4 * SQR(BOA_ij - 1.5) )
* EXP( -p_coa4 * SQR(BOA_jk - 1.5) )
* EXP( -p_coa3 * SQR(workspace->total_bond_order[i] - BOA_ij) )
* EXP( -p_coa3 * SQR(workspace->total_bond_order[k] - BOA_jk) )
/ (1.0 + exp_coa2);
data->my_en.e_coa += e_coa;
CEcoa1 = -2.0 * p_coa4 * (BOA_ij - 1.5) * e_coa;
CEcoa2 = -2.0 * p_coa4 * (BOA_jk - 1.5) * e_coa;
CEcoa3 = -p_coa2 * exp_coa2 * e_coa / (1.0 + exp_coa2);
CEcoa4 = -2.0 * p_coa3 * (workspace->total_bond_order[i] - BOA_ij) * e_coa;
CEcoa5 = -2.0 * p_coa3 * (workspace->total_bond_order[k] - BOA_jk) * e_coa;
/* calculate force contributions */
bo_ij->Cdbo += CEval1 + CEpen2 + (CEcoa1 - CEcoa4);
bo_jk->Cdbo += CEval2 + CEpen3 + (CEcoa2 - CEcoa5);
workspace->CdDelta[j] += ((CEval3 + CEval7) + CEpen1 + CEcoa3);
workspace->CdDelta[i] += CEcoa4;
workspace->CdDelta[k] += CEcoa5;
for ( t = start_j; t < end_j; ++t )
{
pbond_jt = &bond_list->bond_list[t];
bo_jt = &pbond_jt->bo_data;
temp_bo_jt = bo_jt->BO;
temp = CUBE( temp_bo_jt );
pBOjt7 = temp * temp * temp_bo_jt;
bo_jt->Cdbo += CEval6 * pBOjt7;
bo_jt->Cdbopi += CEval5;
bo_jt->Cdbopi2 += CEval5;
}
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 );
}
else
{
/* terms not related to bond order derivatives are
* added directly into forces and pressure vector/tensor */
rvec_Scale( force, CEval8, p_ijk->dcos_di );
rvec_Add( workspace->f[i], force );
rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
rvec_Add( data->my_ext_press, ext_press );
rvec_ScaledAdd( workspace->f[j], CEval8, p_ijk->dcos_dj );
rvec_Scale( force, CEval8, p_ijk->dcos_dk );
rvec_Add( workspace->f[k], force );
rvec_iMultiply( ext_press, pbond_jk->rel_box, force );
rvec_Add( data->my_ext_press, ext_press );
}
pbond_jt = &bond_list->bond_list[t];
bo_jt = &pbond_jt->bo_data;
temp_bo_jt = bo_jt->BO;
temp = CUBE( temp_bo_jt );
pBOjt7 = temp * temp * temp_bo_jt;
bo_jt->Cdbo += CEval6 * pBOjt7;
bo_jt->Cdbopi += CEval5;
bo_jt->Cdbopi2 += CEval5;
}
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 );
}
else
{
/* terms not related to bond order derivatives are
* added directly into forces and pressure vector/tensor */
rvec_Scale( force, CEval8, p_ijk->dcos_di );
rvec_Add( workspace->f[i], force );
rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
rvec_Add( data->my_ext_press, ext_press );
rvec_ScaledAdd( workspace->f[j], CEval8, p_ijk->dcos_dj );
rvec_Scale( force, CEval8, p_ijk->dcos_dk );
rvec_Add( workspace->f[k], force );
rvec_iMultiply( ext_press, pbond_jk->rel_box, force );
rvec_Add( data->my_ext_press, ext_press );
}
#if defined(TEST_ENERGY)
/*fprintf( out_control->eval, "%12.8f%12.8f%12.8f%12.8f\n",
p_val3, p_val4, BOA_ij, BOA_jk );
fprintf(out_control->eval, "%13.8f%13.8f%13.8f%13.8f%13.8f\n",
workspace->Delta_e[j], workspace->vlpex[j],
dSBO1, dSBO2, vlpadj );
fprintf( out_control->eval, "%12.8f%12.8f%12.8f%12.8f\n",
f7_ij, f7_jk, f8_Dj, expval12theta );
fprintf( out_control->eval,
"%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f\n",
CEval1, CEval2, CEval3, CEval4,
CEval5, CEval6, CEval7, CEval8 );
fprintf( out_control->eval,
"%12.8f%12.8f%12.8f\n%12.8f%12.8f%12.8f\n%12.8f%12.8f%12.8f\n",
p_ijk->dcos_di[0]/sin_theta, p_ijk->dcos_di[1]/sin_theta,
p_ijk->dcos_di[2]/sin_theta,
p_ijk->dcos_dj[0]/sin_theta, p_ijk->dcos_dj[1]/sin_theta,
p_ijk->dcos_dj[2]/sin_theta,
p_ijk->dcos_dk[0]/sin_theta, p_ijk->dcos_dk[1]/sin_theta,
p_ijk->dcos_dk[2]/sin_theta);
fprintf( out_control->eval,
"%6d%6d%6d%15.8f%15.8f\n",
/*fprintf( out_control->eval, "%12.8f%12.8f%12.8f%12.8f\n",
p_val3, p_val4, BOA_ij, BOA_jk );
fprintf(out_control->eval, "%13.8f%13.8f%13.8f%13.8f%13.8f\n",
workspace->Delta_e[j], workspace->vlpex[j],
dSBO1, dSBO2, vlpadj );
fprintf( out_control->eval, "%12.8f%12.8f%12.8f%12.8f\n",
f7_ij, f7_jk, f8_Dj, expval12theta );
fprintf( out_control->eval,
"%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f\n",
CEval1, CEval2, CEval3, CEval4,
CEval5, CEval6, CEval7, CEval8 );
fprintf( out_control->eval,
"%12.8f%12.8f%12.8f\n%12.8f%12.8f%12.8f\n%12.8f%12.8f%12.8f\n",
p_ijk->dcos_di[0]/sin_theta, p_ijk->dcos_di[1]/sin_theta,
p_ijk->dcos_di[2]/sin_theta,
p_ijk->dcos_dj[0]/sin_theta, p_ijk->dcos_dj[1]/sin_theta,
p_ijk->dcos_dj[2]/sin_theta,
p_ijk->dcos_dk[0]/sin_theta, p_ijk->dcos_dk[1]/sin_theta,
p_ijk->dcos_dk[2]/sin_theta);
fprintf( out_control->eval,
"%6d%6d%6d%15.8f%15.8f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), e_ang );*/
fprintf( out_control->eval,
//"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), theta_0, BOA_ij, BOA_jk,
e_ang, data->my_en.e_ang );
fprintf( out_control->epen,
//"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), e_ang );*/
fprintf( out_control->eval,
//"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), theta_0, BOA_ij, BOA_jk,
e_ang, data->my_en.e_ang );
fprintf( out_control->epen,
//"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), BOA_ij, BOA_jk, e_pen,
data->my_en.e_pen );
fprintf( out_control->ecoa,
//"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), BOA_ij, BOA_jk,
e_coa, data->my_en.e_coa );
RAD2DEG(theta), BOA_ij, BOA_jk, e_pen,
data->my_en.e_pen );
fprintf( out_control->ecoa,
//"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), BOA_ij, BOA_jk,
e_coa, data->my_en.e_coa );
#endif
#if defined(TEST_FORCES)
/* angle forces */
Add_dBO( system, lists, j, pi, CEval1, workspace->f_ang );
Add_dBO( system, lists, j, pk, CEval2, workspace->f_ang );
Add_dDelta( system, lists, j,
CEval3 + CEval7, workspace->f_ang );
for ( t = start_j; t < end_j; ++t )
{
pbond_jt = &bond_list->bond_list[t];
bo_jt = &pbond_jt->bo_data;
temp_bo_jt = bo_jt->BO;
temp = CUBE( temp_bo_jt );
pBOjt7 = temp * temp * temp_bo_jt;
Add_dBO( system, lists, j, t, pBOjt7 * CEval6,
workspace->f_ang );
Add_dBOpinpi2( system, lists, j, t, CEval5, CEval5,
workspace->f_ang, workspace->f_ang );
}
rvec_ScaledAdd( workspace->f_ang[i], CEval8, p_ijk->dcos_di );
rvec_ScaledAdd( workspace->f_ang[j], CEval8, p_ijk->dcos_dj );
rvec_ScaledAdd( workspace->f_ang[k], CEval8, p_ijk->dcos_dk );
/* end angle forces */
/* penalty forces */
Add_dDelta( system, lists, j, CEpen1, workspace->f_pen );
Add_dBO( system, lists, j, pi, CEpen2, workspace->f_pen );
Add_dBO( system, lists, j, pk, CEpen3, workspace->f_pen );
/* end penalty forces */
/* coalition forces */
Add_dBO( system, lists, j, pi, CEcoa1 - CEcoa4,
workspace->f_coa );
Add_dBO( system, lists, j, pk, CEcoa2 - CEcoa5,
workspace->f_coa );
Add_dDelta( system, lists, j, CEcoa3, workspace->f_coa );
Add_dDelta( system, lists, i, CEcoa4, workspace->f_coa );
Add_dDelta( system, lists, k, CEcoa5, workspace->f_coa );
/* end coalition forces */
#endif
/* angle forces */
Add_dBO( system, lists, j, pi, CEval1, workspace->f_ang );
Add_dBO( system, lists, j, pk, CEval2, workspace->f_ang );
Add_dDelta( system, lists, j,
CEval3 + CEval7, workspace->f_ang );
for ( t = start_j; t < end_j; ++t )
{
pbond_jt = &bond_list->bond_list[t];
bo_jt = &pbond_jt->bo_data;
temp_bo_jt = bo_jt->BO;
temp = CUBE( temp_bo_jt );
pBOjt7 = temp * temp * temp_bo_jt;
Add_dBO( system, lists, j, t, pBOjt7 * CEval6,
workspace->f_ang );
Add_dBOpinpi2( system, lists, j, t, CEval5, CEval5,
workspace->f_ang, workspace->f_ang );
}
rvec_ScaledAdd( workspace->f_ang[i], CEval8, p_ijk->dcos_di );
rvec_ScaledAdd( workspace->f_ang[j], CEval8, p_ijk->dcos_dj );
rvec_ScaledAdd( workspace->f_ang[k], CEval8, p_ijk->dcos_dk );
/* end angle forces */
/* penalty forces */
Add_dDelta( system, lists, j, CEpen1, workspace->f_pen );
Add_dBO( system, lists, j, pi, CEpen2, workspace->f_pen );
Add_dBO( system, lists, j, pk, CEpen3, workspace->f_pen );
/* end penalty forces */
/* coalition forces */
Add_dBO( system, lists, j, pi, CEcoa1 - CEcoa4,
workspace->f_coa );
Add_dBO( system, lists, j, pk, CEcoa2 - CEcoa5,
workspace->f_coa );
Add_dDelta( system, lists, j, CEcoa3, workspace->f_coa );
Add_dDelta( system, lists, i, CEcoa4, workspace->f_coa );
Add_dDelta( system, lists, k, CEcoa5, workspace->f_coa );
/* end coalition forces */
#endif
}
}
}
......
sPuReMD/src/nonbonded.c
+ 3
− 3
View file @ 20ece267
......@@ -164,7 +164,8 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
dTap = dTap * r_ij + workspace->Tap[1];
/* vdWaals Calculations */
if ( system->reax_param.gp.vdw_type == 1 || system->reax_param.gp.vdw_type == 3 )
if ( system->reax_param.gp.vdw_type == 1
|| system->reax_param.gp.vdw_type == 3 )
{
/* shielding */
powr_vdW1 = POW( r_ij, p_vdW1 );
......@@ -221,8 +222,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
#endif
/* Coulomb Calculations */
dr3gamij_1 = r_ij * r_ij * r_ij
+ POW( twbp->gamma, -3.0 );
dr3gamij_1 = r_ij * r_ij * r_ij + POW( twbp->gamma, -3.0 );
dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
e_clb = C_ELE * (system->atoms[i].q * system->atoms[j].q) / dr3gamij_3;
e_ele = self_coef * (e_clb * Tap);
......
sPuReMD/src/valence_angles.c
+ 7
− 7
View file @ 20ece267
......@@ -64,7 +64,7 @@ void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
sqr_d_ji = SQR( d_ji );
sqr_d_jk = SQR( d_jk );
inv_dists = 1.0 / (d_ji * d_jk);
inv_dists3 = POW( inv_dists, 3 );
inv_dists3 = POW( inv_dists, 3.0 );
dot_dvecs = rvec_Dot( dvec_ji, dvec_jk );
Cdot_inv3 = dot_dvecs * inv_dists3;
......@@ -73,7 +73,7 @@ void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
(*dcos_theta_di)[t] = dvec_jk[t] * inv_dists
- Cdot_inv3 * sqr_d_jk * dvec_ji[t];
(*dcos_theta_dj)[t] = -(dvec_jk[t] + dvec_ji[t]) * inv_dists
(*dcos_theta_dj)[t] = -1.0 * (dvec_jk[t] + dvec_ji[t]) * inv_dists
+ Cdot_inv3 * ( sqr_d_jk * dvec_ji[t] + sqr_d_ji * dvec_jk[t] );
(*dcos_theta_dk)[t] = dvec_ji[t] * inv_dists
......@@ -377,13 +377,13 @@ void Valence_Angles( reax_system *system, control_params *control,
f8_Dj = p_val5 - (p_val5 - 1.0) * (2.0 + expval6) / trm8;
Cf8j = ( (1.0 - p_val5) / SQR(trm8) )
* (p_val6 * expval6 * trm8
- (2.0 + expval6) * ( p_val6 * expval6 - p_val7 * expval7 ));
- (2.0 + expval6) * ( p_val6 * expval6 - p_val7 * expval7 ));
theta_0 = 180.0 - theta_00 * (1.0 - EXP(-p_val10 * (2.0 - SBO2)));
theta_0 = DEG2RAD( theta_0 );
expval2theta = p_val1 * EXP(-p_val2 * SQR(theta_0 - theta));
if ( p_val1 >= 0 )
if ( p_val1 >= 0.0 )
{
expval12theta = p_val1 - expval2theta;
}
......@@ -477,15 +477,15 @@ void Valence_Angles( reax_system *system, control_params *control,
#if defined(_OPENMP)
// #pragma omp atomic
#endif
bo_ij->Cdbo += (CEval1 + CEpen2 + (CEcoa1 - CEcoa4));
bo_ij->Cdbo += CEval1 + CEpen2 + (CEcoa1 - CEcoa4);
#if defined(_OPENMP)
// #pragma omp atomic
#endif
bo_jk->Cdbo += (CEval2 + CEpen3 + (CEcoa2 - CEcoa5));
bo_jk->Cdbo += CEval2 + CEpen3 + (CEcoa2 - CEcoa5);
#if defined(_OPENMP)
// #pragma omp atomic
#endif
workspace->CdDelta[j] += ((CEval3 + CEval7) + CEpen1 + CEcoa3);
workspace->CdDelta[j] += (CEval3 + CEval7) + CEpen1 + CEcoa3;
#if defined(_OPENMP)
// #pragma omp atomic
#endif
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment