Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78284811
reaxc_valence_angles.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Aug 19, 14:55
Size
17 KB
Mime Type
text/x-c
Expires
Wed, Aug 21, 14:55 (2 d)
Engine
blob
Format
Raw Data
Handle
20021043
Attached To
rLAMMPS lammps
reaxc_valence_angles.cpp
View Options
/*----------------------------------------------------------------------
PuReMD - Purdue ReaxFF Molecular Dynamics Program
Copyright (2010) Purdue University
Hasan Metin Aktulga, haktulga@cs.purdue.edu
Joseph Fogarty, jcfogart@mail.usf.edu
Sagar Pandit, pandit@usf.edu
Ananth Y Grama, ayg@cs.purdue.edu
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details:
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#if defined(PURE_REAX)
#include "valence_angles.h"
#include "bond_orders.h"
#include "list.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_valence_angles.h"
#include "reaxc_bond_orders.h"
#include "reaxc_list.h"
#include "reaxc_vector.h"
#endif
/* calculates the theta angle between i-j-k */
void Calculate_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
real *theta, real *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;
(*theta) = acos( *cos_theta );
}
/* calculates the derivative of the cosine of the angle between i-j-k */
void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
rvec* dcos_theta_di,
rvec* dcos_theta_dj,
rvec* dcos_theta_dk )
{
int t;
real sqr_d_ji = SQR(d_ji);
real sqr_d_jk = SQR(d_jk);
real inv_dists = 1.0 / (d_ji * d_jk);
real inv_dists3 = pow( inv_dists, 3 );
real dot_dvecs = Dot( dvec_ji, dvec_jk, 3 );
real Cdot_inv3 = dot_dvecs * inv_dists3;
for( t = 0; t < 3; ++t ) {
(*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 +
Cdot_inv3 * ( sqr_d_jk * dvec_ji[t] + sqr_d_ji * dvec_jk[t] );
(*dcos_theta_dk)[t] = dvec_ji[t] * inv_dists -
Cdot_inv3 * sqr_d_ji * dvec_jk[t];
}
}
/* this is a 3-body interaction in which the main role is
played by j which sits in the middle of the other two. */
void Valence_Angles( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control )
{
int i, j, pi, k, pk, t;
int type_i, type_j, type_k;
int start_j, end_j, start_pk, end_pk;
int cnt, num_thb_intrs;
real temp, temp_bo_jt, pBOjt7;
real p_val1, p_val2, p_val3, p_val4, p_val5;
real p_val6, p_val7, p_val8, p_val9, p_val10;
real p_pen1, p_pen2, p_pen3, p_pen4;
real p_coa1, p_coa2, p_coa3, p_coa4;
real trm8, expval6, expval7, expval2theta, expval12theta, exp3ij, exp3jk;
real exp_pen2ij, exp_pen2jk, exp_pen3, exp_pen4, trm_pen34, exp_coa2;
real dSBO1, dSBO2, SBO, SBO2, CSBO2, SBOp, prod_SBO, vlpadj;
real CEval1, CEval2, CEval3, CEval4, CEval5, CEval6, CEval7, CEval8;
real CEpen1, CEpen2, CEpen3;
real e_ang, e_coa, e_pen;
real CEcoa1, CEcoa2, CEcoa3, CEcoa4, CEcoa5;
real Cf7ij, Cf7jk, Cf8j, Cf9j;
real f7_ij, f7_jk, f8_Dj, f9_Dj;
real Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta;
real r_ij, r_jk;
real BOA_ij, BOA_jk;
rvec force, ext_press;
// rtensor temp_rtensor, total_rtensor;
three_body_header *thbh;
three_body_parameters *thbp;
three_body_interaction_data *p_ijk, *p_kji;
bond_data *pbond_ij, *pbond_jk, *pbond_jt;
bond_order_data *bo_ij, *bo_jk, *bo_jt;
reax_list *bonds = (*lists) + BONDS;
reax_list *thb_intrs = (*lists) + THREE_BODIES;
/* global parameters used in these calculations */
p_val6 = system->reax_param.gp.l[14];
p_val8 = system->reax_param.gp.l[33];
p_val9 = system->reax_param.gp.l[16];
p_val10 = system->reax_param.gp.l[17];
num_thb_intrs = 0;
for( j = 0; j < system->N; ++j ) {
// fprintf( out_control->eval, "j: %d\n", j );
type_j = system->my_atoms[j].type;
start_j = Start_Index(j, bonds);
end_j = End_Index(j, bonds);
p_val3 = system->reax_param.sbp[ type_j ].p_val3;
p_val5 = system->reax_param.sbp[ type_j ].p_val5;
SBOp = 0, prod_SBO = 1;
for( t = start_j; t < end_j; ++t ) {
bo_jt = &(bonds->select.bond_list[t].bo_data);
SBOp += (bo_jt->BO_pi + bo_jt->BO_pi2);
temp = SQR( bo_jt->BO );
temp *= temp;
temp *= temp;
prod_SBO *= exp( -temp );
}
/* modifications to match Adri's code - 09/01/09 */
if( workspace->vlpex[j] >= 0 ){
vlpadj = 0;
dSBO2 = prod_SBO - 1;
}
else{
vlpadj = workspace->nlp[j];
dSBO2 = (prod_SBO - 1) * (1 - p_val8 * workspace->dDelta_lp[j]);
}
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 )
SBO2 = 0, CSBO2 = 0;
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 ) {
SBO2 = 2 - pow( 2-SBO, p_val9 );
CSBO2 = p_val9 * pow( 2 - SBO, p_val9 - 1 );
}
else
SBO2 = 2, CSBO2 = 0;
expval6 = exp( p_val6 * workspace->Delta_boc[j] );
for( pi = start_j; pi < end_j; ++pi ) {
Set_Start_Index( pi, num_thb_intrs, thb_intrs );
pbond_ij = &(bonds->select.bond_list[pi]);
bo_ij = &(pbond_ij->bo_data);
BOA_ij = bo_ij->BO - control->thb_cut;
if( BOA_ij/*bo_ij->BO*/ > 0.0 &&
( j < system->n || pbond_ij->nbr < system->n ) ) {
i = pbond_ij->nbr;
r_ij = pbond_ij->d;
type_i = system->my_atoms[i].type;
// fprintf( out_control->eval, "i: %d\n", i );
/* first copy 3-body intrs from previously computed ones where i>k.
in the second for-loop below,
we compute only new 3-body intrs where i < k */
for( pk = start_j; pk < pi; ++pk ) {
// fprintf( out_control->eval, "pk: %d\n", pk );
start_pk = Start_Index( pk, thb_intrs );
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 ) {
p_ijk = &(thb_intrs->select.three_body_list[num_thb_intrs] );
p_kji = &(thb_intrs->select.three_body_list[t]);
p_ijk->thb = bonds->select.bond_list[pk].nbr;
p_ijk->pthb = pk;
p_ijk->theta = p_kji->theta;
rvec_Copy( p_ijk->dcos_di, p_kji->dcos_dk );
rvec_Copy( p_ijk->dcos_dj, p_kji->dcos_dj );
rvec_Copy( p_ijk->dcos_dk, p_kji->dcos_di );
++num_thb_intrs;
break;
}
}
/* and this is the second for loop mentioned above */
for( pk = pi+1; pk < end_j; ++pk ) {
pbond_jk = &(bonds->select.bond_list[pk]);
bo_jk = &(pbond_jk->bo_data);
BOA_jk = bo_jk->BO - control->thb_cut;
k = pbond_jk->nbr;
type_k = system->my_atoms[k].type;
p_ijk = &( thb_intrs->select.three_body_list[num_thb_intrs] );
Calculate_Theta( pbond_ij->dvec, pbond_ij->d,
pbond_jk->dvec, pbond_jk->d,
&theta, &cos_theta );
Calculate_dCos_Theta( pbond_ij->dvec, pbond_ij->d,
pbond_jk->dvec, pbond_jk->d,
&(p_ijk->dcos_di), &(p_ijk->dcos_dj),
&(p_ijk->dcos_dk) );
p_ijk->thb = k;
p_ijk->pthb = pk;
p_ijk->theta = theta;
sin_theta = sin( theta );
if( sin_theta < 1.0e-5 )
sin_theta = 1.0e-5;
++num_thb_intrs;
if( (j < system->n) && (BOA_jk > 0.0) &&
(bo_ij->BO * bo_jk->BO > SQR(control->thb_cut)/*0*/) ) {
r_jk = pbond_jk->d;
thbh = &( system->reax_param.thbp[ type_i ][ type_j ][ type_k ] );
/* if( system->my_atoms[i].orig_id < system->my_atoms[k].orig_id )
fprintf( fval, "%6d %6d %6d %7.3f %7.3f %7.3f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
bo_ij->BO, bo_jk->BO, p_ijk->theta );
else
fprintf( fval, "%6d %6d %6d %7.3f %7.3f %7.3f\n",
system->my_atoms[k].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[i].orig_id,
bo_jk->BO, bo_ij->BO, p_ijk->theta ); */
for( cnt = 0; cnt < thbh->cnt; ++cnt ) {
// fprintf( out_control->eval, "%6d%6d%6d -- exists in thbp\n",
// i+1, j+1, k+1 );
if( fabs(thbh->prm[cnt].p_val1) > 0.001 ) {
thbp = &( thbh->prm[cnt] );
/* 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 = exp( -p_val2 * SQR(theta_0 - theta) );
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;
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_val1 * 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;
data->my_en.e_ang += e_ang =
f7_ij * f7_jk * f8_Dj * expval12theta;
/* END ANGLE ENERGY*/
/* PENALTY ENERGY */
p_pen1 = thbp->p_pen1;
p_pen2 = system->reax_param.gp.l[19];
p_pen3 = system->reax_param.gp.l[20];
p_pen4 = system->reax_param.gp.l[21];
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 );
data->my_en.e_pen += e_pen =
p_pen1 * f9_Dj * exp_pen2ij * exp_pen2jk;
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);
/* END PENALTY ENERGY */
/* COALITION ENERGY */
p_coa1 = thbp->p_coa1;
p_coa2 = system->reax_param.gp.l[2];
p_coa3 = system->reax_param.gp.l[38];
p_coa4 = system->reax_param.gp.l[30];
exp_coa2 = exp( p_coa2 * workspace->Delta_val[j] );
data->my_en.e_coa += e_coa =
p_coa1 / (1. + exp_coa2) *
exp( -p_coa3 * SQR(workspace->total_bond_order[i]-BOA_ij) ) *
exp( -p_coa3 * SQR(workspace->total_bond_order[k]-BOA_jk) ) *
exp( -p_coa4 * SQR(BOA_ij - 1.5) ) *
exp( -p_coa4 * SQR(BOA_jk - 1.5) );
CEcoa1 = -2 * p_coa4 * (BOA_ij - 1.5) * e_coa;
CEcoa2 = -2 * p_coa4 * (BOA_jk - 1.5) * e_coa;
CEcoa3 = -p_coa2 * exp_coa2 * e_coa / (1 + exp_coa2);
CEcoa4 = -2 * p_coa3 *
(workspace->total_bond_order[i]-BOA_ij) * e_coa;
CEcoa5 = -2 * p_coa3 *
(workspace->total_bond_order[k]-BOA_jk) * e_coa;
/* END COALITION ENERGY */
/* FORCES */
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 = &( bonds->select.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;
// fprintf( out_control->eval, "%6d%12.8f\n",
// workspace->reverse_map[bonds->select.bond_list[t].nbr],
// (CEval6 * pBOjt7) );
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 );
}
#ifdef 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",
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 );
#endif
#ifdef 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 = &( bonds->select.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
}
}
}
}
}
Set_End_Index(pi, num_thb_intrs, thb_intrs );
}
}
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 ) {
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 );
}
}
//fprintf( stderr,"%d: Number of angle interactions: %d\n",
// data->step, num_thb_intrs );
#if defined(DEBUG)
fprintf( stderr, "Number of angle interactions: %d\n", num_thb_intrs );
fprintf( stderr,
"Angle Energy: %g\t Penalty Energy: %g\t Coalition Energy: %g\t\n",
data->my_en.e_ang, data->my_en.e_pen, data->my_en.e_coa );
fprintf( stderr, "3body: ext_press (%12.6f %12.6f %12.6f)\n",
data->ext_press[0], data->ext_press[1], data->ext_press[2] );
#endif
}
Event Timeline
Log In to Comment