jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7031 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2011-10-06 14:25:53 +00:00
parent 77371a142f
commit 2e32772d1d
28 changed files with 279 additions and 65 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/USER-REAXC/pair_reax_c.cpp
View File

@ -86,6 +86,7 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp)
system->bndry_cuts.ghost_bond = 0;
system->bndry_cuts.ghost_cutoff = 0;
system->my_atoms = NULL;
system->pair_ptr = this;
fix_reax = NULL;
@ -265,7 +266,7 @@ void PairReaxC::coeff( int nargs, char **args )
void PairReaxC::init_style( )
{
if (!atom->q_flag) error->all(FLERR,"Pair reax/c requires atom attribute q");
firstwarn = 1;
// firstwarn = 1;
int iqeq;
for (iqeq = 0; iqeq < modify->nfix; iqeq++)
@ -395,12 +396,12 @@ void PairReaxC::compute(int eflag, int vflag)
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global = vflag_global = 0;
if ((eflag_atom || vflag_atom) && firstwarn) {
/* if ((eflag_atom || vflag_atom) && firstwarn) {
firstwarn = 0;
if (comm->me == 0)
error->warning(FLERR,"Pair reax/c cannot yet compute "
"per-atom energy or stress");
}
} */
if (vflag_global) control->virial = 1;
else control->virial = 0;
@ -455,8 +456,8 @@ void PairReaxC::compute(int eflag, int vflag)
ecoul += data->my_en.e_ele;
ecoul += data->my_en.e_pol;
eng_vdwl += evdwl;
eng_coul += ecoul;
// eng_vdwl += evdwl;
// eng_coul += ecoul;
// Store the different parts of the energy
// in a list for output by compute pair command

4
src/USER-REAXC/reaxc_allocate.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "allocate.h"
#include "list.h"
@ -885,7 +885,7 @@ void ReAllocate( reax_system *system, control_params *control,
if( system->numH >= DANGER_ZONE * system->Hcap ||
(0 && system->numH <= LOOSE_ZONE * system->Hcap) ) {
Hflag = 1;
system->Hcap = MAX( system->numH * SAFER_ZONE, MIN_CAP );
system->Hcap = int(MAX( system->numH * SAFER_ZONE, MIN_CAP ));
}
if( Hflag || realloc->hbonds ) {

2
src/USER-REAXC/reaxc_basic_comm.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "basic_comm.h"
#include "vector.h"

48
src/USER-REAXC/reaxc_bond_orders.cpp
View File

@ -24,6 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#include "reaxc_types.h"
#if defined(PURE_REAX)
#include "bond_orders.h"
@ -544,7 +545,7 @@ void Add_dBond_to_Forces_NPT( int i, int pj, simulation_data *data,
void Add_dBond_to_Forces( int i, int pj,
void Add_dBond_to_Forces( reax_system *system, int i, int pj,
storage *workspace, reax_list **lists )
{
reax_list *bonds = (*lists) + BONDS;
@ -553,6 +554,10 @@ void Add_dBond_to_Forces( int i, int pj,
dbond_coefficients coef;
int pk, k, j;
/* Virial Tallying variables */
int ii;
real f_scaler, fi_tmp[3], fj_tmp[3], fk_tmp[3];
/* Initializations */
nbr_j = &(bonds->select.bond_list[pj]);
j = nbr_j->nbr;
@ -589,6 +594,13 @@ void Add_dBond_to_Forces( int i, int pj,
rvec_ScaledAdd( workspace->f[k], -coef.C3dbopi, nbr_k->bo_data.dBOp );
/*3rd, dBOpi2*/
rvec_ScaledAdd( workspace->f[k], -coef.C3dbopi2, nbr_k->bo_data.dBOp );
// tally into per-atom virial
if( system->vflag_atom ) {
f_scaler = -(coef.C2dbo + coef.C2dDelta + coef.C3dbopi + coef.C3dbopi2);
rvec_Scale( fk_tmp, -f_scaler, nbr_k->bo_data.dBOp);
system->pair_ptr->v_tally(k, fk_tmp);
}
}
/*1st, dBO*/
@ -628,6 +640,13 @@ void Add_dBond_to_Forces( int i, int pj,
rvec_ScaledAdd( workspace->f[k], -coef.C4dbopi, nbr_k->bo_data.dBOp );
/*4th, dBOpi2*/
rvec_ScaledAdd( workspace->f[k], -coef.C4dbopi2, nbr_k->bo_data.dBOp );
// tally into per-atom virial
if( system->vflag_atom ) {
f_scaler = -(coef.C3dbo + coef.C3dDelta + coef.C4dbopi + coef.C4dbopi2);
rvec_Scale( fk_tmp, -f_scaler, nbr_k->bo_data.dBOp);
system->pair_ptr->v_tally(k, fk_tmp);
}
}
/*1st,dBO*/
@ -653,6 +672,33 @@ void Add_dBond_to_Forces( int i, int pj,
rvec_ScaledAdd( workspace->f[j], -coef.C2dbopi2, bo_ij->dBOp );
/*3rd, dBOpi2*/
rvec_ScaledAdd( workspace->f[j], coef.C4dbopi2, workspace->dDeltap_self[j] );
if( system->vflag_atom) {
// forces on i
f_scaler = coef.C1dbo + coef.C1dDelta + coef.C2dbopi + coef.C2dbopi2;
rvec_Scale( fi_tmp, -f_scaler, bo_ij->dBOp);
f_scaler = coef.C2dbo + coef.C2dDelta + coef.C3dbopi + coef.C3dbopi2;
rvec_ScaledAdd( fi_tmp, -f_scaler, workspace->dDeltap_self[i]);
f_scaler = coef.C1dbopi ;
rvec_ScaledAdd( fi_tmp, -f_scaler, bo_ij->dln_BOp_pi);
f_scaler = coef.C1dbopi2 ;
rvec_ScaledAdd( fi_tmp, -f_scaler, bo_ij->dln_BOp_pi2);
// forces on j
f_scaler = -(coef.C1dbo + coef.C1dDelta + coef.C2dbopi + coef.C2dbopi2);
rvec_Scale( fj_tmp, -f_scaler, bo_ij->dBOp);
f_scaler = coef.C3dbo + coef.C3dDelta + coef.C4dbopi + coef.C4dbopi2;
rvec_ScaledAdd( fj_tmp, -f_scaler, workspace->dDeltap_self[j]);
f_scaler = -coef.C1dbopi ;
rvec_ScaledAdd( fj_tmp, -f_scaler, bo_ij->dln_BOp_pi);
f_scaler = -coef.C1dbopi2 ;
rvec_ScaledAdd( fj_tmp, -f_scaler, bo_ij->dln_BOp_pi2);
system->pair_ptr->v_tally(i, fi_tmp);
system->pair_ptr->v_tally(j, fj_tmp);
}
}

2
src/USER-REAXC/reaxc_bond_orders.h
View File

@ -53,7 +53,7 @@ void Add_dDelta( reax_system*, reax_list**, int, real, rvec* );
void Add_dDelta_to_Forces( reax_system *, reax_list**, int, real );
#endif
void Add_dBond_to_Forces( int, int, storage*, reax_list** );
void Add_dBond_to_Forces( reax_system*, int, int, storage*, reax_list** );
void Add_dBond_to_Forces_NPT( int, int, simulation_data*,
storage*, reax_list** );
int BOp(storage*, reax_list*, real, int, int, far_neighbor_data*,

14
src/USER-REAXC/reaxc_bonds.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "bonds.h"
#include "bond_orders.h"
@ -79,7 +79,7 @@ void Bonds( reax_system *system, control_params *control,
sbp_j = &( system->reax_param.sbp[type_j] );
twbp = &( system->reax_param.tbp[type_i][type_j] );
bo_ij = &( bonds->select.bond_list[pj].bo_data );
/* calculate the constants */
pow_BOs_be2 = pow( bo_ij->BO_s, twbp->p_be2 );
exp_be12 = exp( twbp->p_be1 * ( 1.0 - pow_BOs_be2 ) );
@ -91,7 +91,11 @@ void Bonds( reax_system *system, control_params *control,
-twbp->De_s * bo_ij->BO_s * exp_be12
-twbp->De_p * bo_ij->BO_pi
-twbp->De_pp * bo_ij->BO_pi2;
/* tally into per-atom energy */
if( system->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 */
bo_ij->Cdbo += CEbo;
bo_ij->Cdbopi -= (CEbo + twbp->De_p);
@ -131,6 +135,10 @@ void Bonds( reax_system *system, control_params *control,
decobdboub = -gp10 * exphu * hulpov *
(gp3*exphub1 + 25.0*gp4*exphuov*hulpov*(exphua1+exphub1));
/* tally into per-atom energy */
if( system->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;
workspace->CdDelta[i] += decobdboua;
workspace->CdDelta[j] += decobdboub;

2
src/USER-REAXC/reaxc_control.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "control.h"
#include "tool_box.h"

4
src/USER-REAXC/reaxc_defs.h
View File

@ -40,8 +40,8 @@
#define CUBE(x) ((x)*(x)*(x))
#define DEG2RAD(a) ((a)*constPI/180.0)
#define RAD2DEG(a) ((a)*180.0/constPI)
#define MAX(x,y) (((x) > (y)) ? (x) : (y))
#define MIN(x,y) (((x) < (y)) ? (x) : (y))
// #define MAX(x,y) (((x) > (y)) ? (x) : (y))
// #define MIN(x,y) (((x) < (y)) ? (x) : (y))
#define MAX3(x,y,z) MAX( MAX(x,y), z)
#define constPI 3.14159265

2
src/USER-REAXC/reaxc_ffield.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "ffield.h"
#include "tool_box.h"

4
src/USER-REAXC/reaxc_forces.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "forces.h"
#include "bond_orders.h"
@ -152,7 +152,7 @@ void Compute_Total_Force( reax_system *system, control_params *control,
for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
if( i < bonds->select.bond_list[pj].nbr ) {
if( control->virial == 0 )
Add_dBond_to_Forces( i, pj, workspace, lists );
Add_dBond_to_Forces( system, i, pj, workspace, lists );
else
Add_dBond_to_Forces_NPT( i, pj, data, workspace, lists );
}

25
src/USER-REAXC/reaxc_hydrogen_bonds.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "hydrogen_bonds.h"
#include "bond_orders.h"
@ -64,6 +64,9 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
bond_data *bond_list;
hbond_data *hbond_list;
// tally variables
real fi_tmp[3], fk_tmp[3], delij[3], delkj[3];
bonds = (*lists) + BONDS;
bond_list = bonds->select.bond_list;
hbonds = (*lists) + HBONDS;
@ -134,14 +137,16 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
exp_hb2 = exp( -hbp->p_hb2 * bo_ij->BO );
exp_hb3 = exp( -hbp->p_hb3 * ( hbp->r0_hb / r_jk +
r_jk / hbp->r0_hb - 2.0 ) );
data->my_en.e_hb += e_hb =
hbp->p_hb1 * (1.0 - exp_hb2) * exp_hb3 * sin_xhz4;
CEhb1 = hbp->p_hb1 * hbp->p_hb2 * exp_hb2 * exp_hb3 * sin_xhz4;
CEhb2 = -hbp->p_hb1/2.0 * (1.0 - exp_hb2) * exp_hb3 * cos_xhz1;
CEhb3 = -hbp->p_hb3 *
(-hbp->r0_hb / SQR(r_jk) + 1.0 / hbp->r0_hb) * e_hb;
fprintf(stderr," H bond called \n");
/*fprintf( stdout,
"%6d%6d%6d%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n",
@ -186,6 +191,20 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
rvec_ScaledAdd( data->my_ext_press, 1.0, ext_press );
}
/* tally into per-atom virials */
if (system->vflag_atom || system->evflag) {
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
-1., system->my_atoms[j].x );
rvec_ScaledSum( delkj, 1., system->my_atoms[k].x,
-1., system->my_atoms[j].x );
rvec_Scale(fi_tmp, CEhb2, dcos_theta_di);
rvec_Scale(fk_tmp, CEhb2, dcos_theta_dk);
rvec_ScaledAdd(fk_tmp, CEhb3/r_jk, dvec_jk);
system->pair_ptr->ev_tally3(i,j,k,e_hb,0.0,fi_tmp,fk_tmp,delij,delkj);
}
#ifdef TEST_ENERGY
/* fprintf( out_control->ehb,
"%24.15e%24.15e%24.15e\n%24.15e%24.15e%24.15e\n%24.15e%24.15e%24.15e\n",

2
src/USER-REAXC/reaxc_init_md.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "init_md.h"
#include "allocate.h"

2
src/USER-REAXC/reaxc_io_tools.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "io_tools.h"
#include "basic_comm.h"

2
src/USER-REAXC/reaxc_list.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "list.h"
#include "tool_box.h"

2
src/USER-REAXC/reaxc_lookup.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "lookup.h"
#include "nonbonded.h"

24
src/USER-REAXC/reaxc_multi_body.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "multi_body.h"
#include "bond_orders.h"
@ -53,6 +53,7 @@ void Atom_Energy( reax_system *system, control_params *control,
real e_un, CEunder1, CEunder2, CEunder3, CEunder4;
real p_lp1, p_lp2, p_lp3;
real p_ovun2, p_ovun3, p_ovun4, p_ovun5, p_ovun6, p_ovun7, p_ovun8;
real eng_tmp, f_tmp;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
@ -82,13 +83,17 @@ void Atom_Energy( reax_system *system, control_params *control,
/* calculate the energy */
data->my_en.e_lp += e_lp =
p_lp2 * workspace->Delta_lp[i] * inv_expvd2;
dElp = p_lp2 * inv_expvd2 +
75 * p_lp2 * workspace->Delta_lp[i] * expvd2 * SQR(inv_expvd2);
CElp = dElp * workspace->dDelta_lp[i];
workspace->CdDelta[i] += CElp; // lp - 1st term
/* tally into per-atom energy */
if( system->evflag)
system->pair_ptr->ev_tally(i,i,system->n,1,e_lp,0.0,0.0,0.0,0.0,0.0);
#ifdef TEST_ENERGY
// fprintf( out_control->elp, "%24.15e%24.15e%24.15e%24.15e\n",
// p_lp2, workspace->Delta_lp_temp[i], expvd2, dElp );
@ -115,12 +120,17 @@ void Atom_Energy( reax_system *system, control_params *control,
if( vov3 > 3. ) {
data->my_en.e_lp += e_lph = p_lp3 * SQR(vov3-3.0);
deahu2dbo = 2.*p_lp3*(vov3 - 3.);
deahu2dsbo = 2.*p_lp3*(vov3 - 3.)*(-1. - 0.16*pow(Di, 3.));
bo_ij->Cdbo += deahu2dbo;
workspace->CdDelta[i] += deahu2dsbo;
/* tally into per-atom energy */
if( system->evflag)
system->pair_ptr->ev_tally(i,j,system->n,1,e_lph,0.0,0.0,0.0,0.0,0.0);
#ifdef TEST_ENERGY
fprintf(out_control->elp,"C2cor%6d%6d%12.6f%12.6f%12.6f\n",
system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
@ -207,7 +217,13 @@ void Atom_Energy( reax_system *system, control_params *control,
CEunder3 = CEunder1 * (1.0 - dfvl*workspace->dDelta_lp[i]*inv_exp_ovun1);
CEunder4 = CEunder1 * (dfvl*workspace->Delta_lp_temp[i]) *
p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1) + CEunder2;
/* tally into per-atom energy */
if( system->evflag) {
eng_tmp = e_ov + e_un;
f_tmp = CEover3 + CEunder3;
system->pair_ptr->ev_tally(i,i,system->n,1,eng_tmp,0.0,0.0,0.0,0.0,0.0);
}
/* forces */
workspace->CdDelta[i] += CEover3; // OvCoor - 2nd term

65
src/USER-REAXC/reaxc_nonbonded.cpp
View File

@ -24,6 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#include "reaxc_types.h"
#if defined(PURE_REAX)
#include "nonbonded.h"
@ -37,7 +38,6 @@
#include "reaxc_vector.h"
#endif
void vdW_Coulomb_Energy( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control )
@ -56,6 +56,9 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
reax_list *far_nbrs;
// rtensor temp_rtensor, total_rtensor;
// Tallying variables:
real pe_vdw, f_tmp, delij[3];
natoms = system->n;
far_nbrs = (*lists) + FAR_NBRS;
p_vdW1 = system->reax_param.gp.l[28];
@ -93,6 +96,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
if (flag) {
#endif
r_ij = nbr_pj->d;
twbp = &(system->reax_param.tbp[ system->my_atoms[i].type ]
[ system->my_atoms[j].type ]);
@ -161,13 +165,23 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
data->my_en.e_ele += e_ele =
C_ele * system->my_atoms[i].q * system->my_atoms[j].q * tmp;
CEclmb = C_ele * system->my_atoms[i].q * system->my_atoms[j].q *
( dTap - Tap * r_ij / dr3gamij_1 ) / dr3gamij_3;
// fprintf( fout, "%5d %5d %10.6f %10.6f\n",
// MIN( system->my_atoms[i].orig_id, system->my_atoms[j].orig_id ),
// MAX( system->my_atoms[i].orig_id, system->my_atoms[j].orig_id ),
// CEvd, CEclmb );
/* tally into per-atom energy */
if( system->evflag || system->vflag_atom) {
pe_vdw = Tap * (e_vdW + e_core);
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
-1., system->my_atoms[j].x );
f_tmp = -(CEvd + CEclmb);
system->pair_ptr->ev_tally(i,j,natoms,1,pe_vdw,e_ele,
f_tmp,delij[0],delij[1],delij[2]);
}
/* fprintf(stderr, "%5d %5d %10.6f %10.6f\n",
MIN( system->my_atoms[i].orig_id, system->my_atoms[j].orig_id ),
MAX( system->my_atoms[i].orig_id, system->my_atoms[j].orig_id ),
CEvd, CEclmb ); */
if( control->virial == 0 ) {
rvec_ScaledAdd( workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec );
@ -184,13 +198,13 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
rvec_iMultiply( ext_press, nbr_pj->rel_box, temp );
rvec_Add( data->my_ext_press, ext_press );
// fprintf( stderr, "nonbonded(%d,%d): rel_box (%f %f %f)
// force(%f %f %f) ext_press (%12.6f %12.6f %12.6f)\n",
// i, j, nbr_pj->rel_box[0], nbr_pj->rel_box[1], nbr_pj->rel_box[2],
// temp[0], temp[1], temp[2],
// data->ext_press[0], data->ext_press[1], data->ext_press[2] );
/* fprintf( stderr, "nonbonded(%d,%d): rel_box (%f %f %f)
force(%f %f %f) ext_press (%12.6f %12.6f %12.6f)\n",
i, j, nbr_pj->rel_box[0], nbr_pj->rel_box[1], nbr_pj->rel_box[2],
temp[0], temp[1], temp[2],
data->ext_press[0], data->ext_press[1], data->ext_press[2] ); */
}
#ifdef TEST_ENERGY
// fprintf( out_control->evdw,
// "%12.9f%12.9f%12.9f%12.9f%12.9f%12.9f%12.9f%12.9f\n",
@ -239,6 +253,8 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
real r_ij, base, dif;
real e_vdW, e_ele;
real CEvd, CEclmb, SMALL = 0.0001;
real f_tmp, delij[3];
rvec temp, ext_press;
far_neighbor_data *nbr_pj;
reax_list *far_nbrs;
@ -283,6 +299,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
#endif
j = nbr_pj->nbr;
type_j = system->my_atoms[j].type;
r_ij = nbr_pj->d;
tmin = MIN( type_i, type_j );
tmax = MAX( type_i, type_j );
@ -314,7 +331,17 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
CEclmb = ((t->CEclmb[r].d*dif+t->CEclmb[r].c)*dif+t->CEclmb[r].b)*dif +
t->CEclmb[r].a;
CEclmb *= system->my_atoms[i].q * system->my_atoms[j].q;
/* tally into per-atom energy */
if( system->evflag || system->vflag_atom) {
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
-1., system->my_atoms[j].x );
f_tmp = -(CEvd + CEclmb);
system->pair_ptr->ev_tally(i,j,natoms,1,e_vdW,e_ele,
f_tmp,delij[0],delij[1],delij[2]);
}
if( control->virial == 0 ) {
rvec_ScaledAdd( workspace->f[i], -(CEvd + CEclmb), nbr_pj->dvec );
rvec_ScaledAdd( workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec );
@ -360,16 +387,20 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
void Compute_Polarization_Energy( reax_system *system, simulation_data *data )
{
int i, type_i;
real q;
real q, en_tmp;
data->my_en.e_pol = 0.0;
for( i = 0; i < system->n; i++ ) {
q = system->my_atoms[i].q;
type_i = system->my_atoms[i].type;
data->my_en.e_pol +=
KCALpMOL_to_EV * (system->reax_param.sbp[type_i].chi * q +
(system->reax_param.sbp[type_i].eta / 2.) * SQR(q));
en_tmp = KCALpMOL_to_EV * (system->reax_param.sbp[type_i].chi * q +
(system->reax_param.sbp[type_i].eta / 2.) * SQR(q));
data->my_en.e_pol += en_tmp;
/* tally into per-atom energy */
if( system->evflag)
system->pair_ptr->ev_tally(i,i,system->n,1,0.0,en_tmp,0.0,0.0,0.0,0.0);
}
}

2
src/USER-REAXC/reaxc_reset_tools.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "reset_tools.h"
#include "list.h"

2
src/USER-REAXC/reaxc_system_props.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "system_props.h"
#include "tool_box.h"

2
src/USER-REAXC/reaxc_tool_box.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "tool_box.h"
#elif defined(LAMMPS_REAX)

39
src/USER-REAXC/reaxc_torsion_angles.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "torsion_angles.h"
#include "bond_orders.h"
@ -201,6 +201,10 @@ void Torsion_Angles( reax_system *system, control_params *control,
// char fname[100];
// FILE *ftor;
// Virial tallying variables
real delil[3], deljl[3], delkl[3];
real eng_tmp, f_scaler, fi_tmp[3], fj_tmp[3], fk_tmp[3];
// sprintf( fname, "tor%d.out", system->my_rank );
// ftor = fopen( fname, "w" );
@ -479,6 +483,39 @@ void Torsion_Angles( reax_system *system, control_params *control,
rvec_Add( data->my_ext_press, ext_press );
}
/* tally into per-atom virials */
if( system->vflag_atom || system->evflag) {
// acquire vectors
rvec_ScaledSum( delil, 1., system->my_atoms[l].x,
-1., system->my_atoms[i].x );
rvec_ScaledSum( deljl, 1., system->my_atoms[l].x,
-1., system->my_atoms[j].x );
rvec_ScaledSum( delkl, 1., system->my_atoms[l].x,
-1., system->my_atoms[k].x );
// dcos_theta_ijk
rvec_Scale( fi_tmp, CEtors7 + CEconj4, p_ijk->dcos_dk );
rvec_Scale( fj_tmp, CEtors7 + CEconj4, p_ijk->dcos_dj );
rvec_Scale( fk_tmp, CEtors7 + CEconj4, p_ijk->dcos_di );
// dcos_theta_jkl
rvec_ScaledAdd( fj_tmp, CEtors8 + CEconj5, p_jkl->dcos_di );
rvec_ScaledAdd( fk_tmp, CEtors8 + CEconj5, p_jkl->dcos_dj );
// dcos_omega
rvec_ScaledAdd( fi_tmp, CEtors9 + CEconj6, dcos_omega_di );
rvec_ScaledAdd( fj_tmp, CEtors9 + CEconj6, dcos_omega_dj );
rvec_ScaledAdd( fk_tmp, CEtors9 + CEconj6, dcos_omega_dk );
// tally
eng_tmp = e_tor + e_con;
if( system->evflag)
system->pair_ptr->ev_tally(j,k,natoms,1,eng_tmp,0.0,0.0,0.0,0.0,0.0);
if( system->vflag_atom)
system->pair_ptr->v_tally4(i,j,k,l,fi_tmp,fj_tmp,fk_tmp,delil,deljl,delkl);
}
#ifdef TEST_ENERGY
/* fprintf( out_control->etor,
"%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f\n",

2
src/USER-REAXC/reaxc_traj.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "traj.h"
#include "list.h"

6
src/USER-REAXC/reaxc_types.h
View File

@ -456,7 +456,7 @@ typedef struct
real ghost_cutoff;
} boundary_cutoff;
using LAMMPS_NS::Pair;
typedef struct
{
@ -474,6 +474,10 @@ typedef struct
boundary_cutoff bndry_cuts;
reax_atom *my_atoms;
int evflag;
int vflag_atom;
class Pair *pair_ptr;
} reax_system;

32
src/USER-REAXC/reaxc_valence_angles.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "valence_angles.h"
#include "bond_orders.h"
@ -106,6 +106,10 @@ void Valence_Angles( reax_system *system, control_params *control,
rvec force, ext_press;
// rtensor temp_rtensor, total_rtensor;
// Tallying variables
real eng_tmp, f_scaler, fi_tmp[3], fj_tmp[3], fk_tmp[3];
real delij[3], delkj[3];
three_body_header *thbh;
three_body_parameters *thbp;
three_body_interaction_data *p_ijk, *p_kji;
@ -114,7 +118,6 @@ void Valence_Angles( reax_system *system, control_params *control,
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];
@ -123,7 +126,7 @@ void Valence_Angles( reax_system *system, control_params *control,
num_thb_intrs = 0;
for( j = 0; j < system->N; ++j ) {
for( j = 0; j < system->N; ++j ) { // Ray: the first one with system->N
// fprintf( out_control->eval, "j: %d\n", j );
type_j = system->my_atoms[j].type;
start_j = Start_Index(j, bonds);
@ -366,6 +369,7 @@ void Valence_Angles( reax_system *system, control_params *control,
CEcoa5 = -2 * p_coa3 *
(workspace->total_bond_order[k]-BOA_jk) * e_coa;
/* END COALITION ENERGY */
/* FORCES */
bo_ij->Cdbo += (CEval1 + CEpen2 + (CEcoa1 - CEcoa4));
@ -390,7 +394,6 @@ void Valence_Angles( reax_system *system, control_params *control,
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 );
@ -411,6 +414,27 @@ void Valence_Angles( reax_system *system, control_params *control,
rvec_iMultiply( ext_press, pbond_jk->rel_box, force );
rvec_Add( data->my_ext_press, ext_press );
}
/* tally into per-atom virials */
if( system->vflag_atom || system->evflag) {
/* Acquire vectors */
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
-1., system->my_atoms[j].x );
rvec_ScaledSum( delkj, 1., system->my_atoms[k].x,
-1., system->my_atoms[j].x );
rvec_Scale( fi_tmp, -CEval8, p_ijk->dcos_di );
rvec_Scale( fj_tmp, -CEval8, p_ijk->dcos_dj );
rvec_Scale( fk_tmp, -CEval8, p_ijk->dcos_dk );
eng_tmp = e_ang + e_pen + e_coa;
if( system->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->vflag_atom)
system->pair_ptr->v_tally3(i,j,k,fi_tmp,fk_tmp,delij,delkj);
}
#ifdef TEST_ENERGY
/*fprintf( out_control->eval, "%12.8f%12.8f%12.8f%12.8f\n",

2
src/USER-REAXC/reaxc_vector.cpp
View File

@ -24,7 +24,7 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxc_types.h"
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "vector.h"
#include "random.h"

1
src/USER-REAXC/reaxc_vector.h
View File

@ -27,6 +27,7 @@
#ifndef __VECTOR_H_
#define __VECTOR_H_
#include "pair.h"
#include "reaxc_types.h"
#include "reaxc_defs.h"

22
src/pair.cpp
View File

@ -759,6 +759,28 @@ void Pair::ev_tally_list(int n, int *list, double ecoul, double *v)
}
}
/* ----------------------------------------------------------------------
tally virial into per-atom accumulators
called by REAX/C potential, newton_pair is always on
fi is magnitude of force on atom i
------------------------------------------------------------------------- */
void Pair::v_tally(int i, double *fi)
{
double v[6];
double **x = atom->x;
v[0] = x[i][0]*fi[0];
v[1] = x[i][1]*fi[1];
v[2] = x[i][2]*fi[2];
v[3] = x[i][0]*fi[1];
v[4] = x[i][0]*fi[1];
v[5] = x[i][1]*fi[2];
vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2];
vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5];
}
/* ----------------------------------------------------------------------
tally virial into per-atom accumulators
called by AIREBO potential, newton_pair is always on

19
src/pair.h
View File

@ -71,6 +71,17 @@ class Pair : protected Pointers {
void init_bitmap(double, double, int, int &, int &, int &, int &);
virtual void modify_params(int, char **);
// need to be public, so can be called by pair_style reaxc
void v_tally(int, double *);
void ev_tally(int, int, int, int, double, double, double,
double, double, double);
void ev_tally3(int, int, int, double, double,
double *, double *, double *, double *);
void v_tally3(int, int, int, double *, double *, double *, double *);
void v_tally4(int, int, int, int, double *, double *, double *,
double *, double *, double *);
// general child-class methods
virtual void compute(int, int) = 0;
@ -125,26 +136,20 @@ class Pair : protected Pointers {
int evflag; // energy,virial settings
int eflag_either,eflag_global,eflag_atom;
int vflag_either,vflag_global,vflag_atom;
int vflag_fdotr;
int maxeatom,maxvatom;
virtual void ev_setup(int, int);
void ev_tally(int, int, int, int, double, double, double,
double, double, double);
void ev_tally_full(int, double, double, double, double, double, double);
void ev_tally_xyz(int, int, int, int, double, double,
double, double, double, double, double, double);
void ev_tally_xyz_full(int, double, double,
double, double, double, double, double, double);
void ev_tally3(int, int, int, double, double,
double *, double *, double *, double *);
void ev_tally4(int, int, int, int, double,
double *, double *, double *, double *, double *, double *);
void ev_tally_list(int, int *, double, double *);
void v_tally2(int, int, double, double *);
void v_tally3(int, int, int, double *, double *, double *, double *);
void v_tally4(int, int, int, int, double *, double *, double *,
double *, double *, double *);
void v_tally_tensor(int, int, int, int,
double, double, double, double, double, double);
void virial_fdotr_compute();