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@13861 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2015-08-10 23:50:55 +00:00
parent 76cfa65725
commit 5f76611ecf
43 changed files with 482 additions and 538 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -777,7 +777,7 @@ double PairReaxC::memory_usage()
// From reaxc_allocate: BO
bytes += 1.0 * system->total_cap * sizeof(reax_atom);
bytes += 19.0 * system->total_cap * sizeof(real);
bytes += 19.0 * system->total_cap * sizeof(double);
bytes += 3.0 * system->total_cap * sizeof(int);
// From reaxc_lists

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

@ -191,14 +191,14 @@ int Allocate_Workspace( reax_system *system, control_params *control,
int i, total_real, total_rvec, local_rvec;
workspace->allocated = 1;
total_real = total_cap * sizeof(real);
total_real = total_cap * sizeof(double);
total_rvec = total_cap * sizeof(rvec);
local_rvec = local_cap * sizeof(rvec);
/* communication storage */
for( i = 0; i < MAX_NBRS; ++i ) {
workspace->tmp_dbl[i] = (real*)
scalloc( total_cap, sizeof(real), "tmp_dbl", comm );
workspace->tmp_dbl[i] = (double*)
scalloc( total_cap, sizeof(double), "tmp_dbl", comm );
workspace->tmp_rvec[i] = (rvec*)
scalloc( total_cap, sizeof(rvec), "tmp_rvec", comm );
workspace->tmp_rvec2[i] = (rvec2*)
@ -208,62 +208,62 @@ int Allocate_Workspace( reax_system *system, control_params *control,
/* bond order related storage */
workspace->within_bond_box = (int*)
scalloc( total_cap, sizeof(int), "skin", comm );
workspace->total_bond_order = (real*) smalloc( total_real, "total_bo", comm );
workspace->Deltap = (real*) smalloc( total_real, "Deltap", comm );
workspace->Deltap_boc = (real*) smalloc( total_real, "Deltap_boc", comm );
workspace->total_bond_order = (double*) smalloc( total_real, "total_bo", comm );
workspace->Deltap = (double*) smalloc( total_real, "Deltap", comm );
workspace->Deltap_boc = (double*) smalloc( total_real, "Deltap_boc", comm );
workspace->dDeltap_self = (rvec*) smalloc( total_rvec, "dDeltap_self", comm );
workspace->Delta = (real*) smalloc( total_real, "Delta", comm );
workspace->Delta_lp = (real*) smalloc( total_real, "Delta_lp", comm );
workspace->Delta_lp_temp = (real*)
workspace->Delta = (double*) smalloc( total_real, "Delta", comm );
workspace->Delta_lp = (double*) smalloc( total_real, "Delta_lp", comm );
workspace->Delta_lp_temp = (double*)
smalloc( total_real, "Delta_lp_temp", comm );
workspace->dDelta_lp = (real*) smalloc( total_real, "dDelta_lp", comm );
workspace->dDelta_lp_temp = (real*)
workspace->dDelta_lp = (double*) smalloc( total_real, "dDelta_lp", comm );
workspace->dDelta_lp_temp = (double*)
smalloc( total_real, "dDelta_lp_temp", comm );
workspace->Delta_e = (real*) smalloc( total_real, "Delta_e", comm );
workspace->Delta_boc = (real*) smalloc( total_real, "Delta_boc", comm );
workspace->Delta_val = (real*) smalloc( total_real, "Delta_val", comm );
workspace->nlp = (real*) smalloc( total_real, "nlp", comm );
workspace->nlp_temp = (real*) smalloc( total_real, "nlp_temp", comm );
workspace->Clp = (real*) smalloc( total_real, "Clp", comm );
workspace->vlpex = (real*) smalloc( total_real, "vlpex", comm );
workspace->Delta_e = (double*) smalloc( total_real, "Delta_e", comm );
workspace->Delta_boc = (double*) smalloc( total_real, "Delta_boc", comm );
workspace->Delta_val = (double*) smalloc( total_real, "Delta_val", comm );
workspace->nlp = (double*) smalloc( total_real, "nlp", comm );
workspace->nlp_temp = (double*) smalloc( total_real, "nlp_temp", comm );
workspace->Clp = (double*) smalloc( total_real, "Clp", comm );
workspace->vlpex = (double*) smalloc( total_real, "vlpex", comm );
workspace->bond_mark = (int*)
scalloc( total_cap, sizeof(int), "bond_mark", comm );
workspace->done_after = (int*)
scalloc( total_cap, sizeof(int), "done_after", comm );
/* QEq storage */
workspace->Hdia_inv = (real*)
scalloc( total_cap, sizeof(real), "Hdia_inv", comm );
workspace->b_s = (real*) scalloc( total_cap, sizeof(real), "b_s", comm );
workspace->b_t = (real*) scalloc( total_cap, sizeof(real), "b_t", comm );
workspace->b_prc = (real*) scalloc( total_cap, sizeof(real), "b_prc", comm );
workspace->b_prm = (real*) scalloc( total_cap, sizeof(real), "b_prm", comm );
workspace->s = (real*) scalloc( total_cap, sizeof(real), "s", comm );
workspace->t = (real*) scalloc( total_cap, sizeof(real), "t", comm );
workspace->droptol = (real*)
scalloc( total_cap, sizeof(real), "droptol", comm );
workspace->Hdia_inv = (double*)
scalloc( total_cap, sizeof(double), "Hdia_inv", comm );
workspace->b_s = (double*) scalloc( total_cap, sizeof(double), "b_s", comm );
workspace->b_t = (double*) scalloc( total_cap, sizeof(double), "b_t", comm );
workspace->b_prc = (double*) scalloc( total_cap, sizeof(double), "b_prc", comm );
workspace->b_prm = (double*) scalloc( total_cap, sizeof(double), "b_prm", comm );
workspace->s = (double*) scalloc( total_cap, sizeof(double), "s", comm );
workspace->t = (double*) scalloc( total_cap, sizeof(double), "t", comm );
workspace->droptol = (double*)
scalloc( total_cap, sizeof(double), "droptol", comm );
workspace->b = (rvec2*) scalloc( total_cap, sizeof(rvec2), "b", comm );
workspace->x = (rvec2*) scalloc( total_cap, sizeof(rvec2), "x", comm );
/* GMRES storage */
workspace->y = (real*) scalloc( RESTART+1, sizeof(real), "y", comm );
workspace->z = (real*) scalloc( RESTART+1, sizeof(real), "z", comm );
workspace->g = (real*) scalloc( RESTART+1, sizeof(real), "g", comm );
workspace->h = (real**) scalloc( RESTART+1, sizeof(real*), "h", comm );
workspace->hs = (real*) scalloc( RESTART+1, sizeof(real), "hs", comm );
workspace->hc = (real*) scalloc( RESTART+1, sizeof(real), "hc", comm );
workspace->v = (real**) scalloc( RESTART+1, sizeof(real*), "v", comm );
workspace->y = (double*) scalloc( RESTART+1, sizeof(double), "y", comm );
workspace->z = (double*) scalloc( RESTART+1, sizeof(double), "z", comm );
workspace->g = (double*) scalloc( RESTART+1, sizeof(double), "g", comm );
workspace->h = (double**) scalloc( RESTART+1, sizeof(double*), "h", comm );
workspace->hs = (double*) scalloc( RESTART+1, sizeof(double), "hs", comm );
workspace->hc = (double*) scalloc( RESTART+1, sizeof(double), "hc", comm );
workspace->v = (double**) scalloc( RESTART+1, sizeof(double*), "v", comm );
for( i = 0; i < RESTART+1; ++i ) {
workspace->h[i] = (real*) scalloc( RESTART+1, sizeof(real), "h[i]", comm );
workspace->v[i] = (real*) scalloc( total_cap, sizeof(real), "v[i]", comm );
workspace->h[i] = (double*) scalloc( RESTART+1, sizeof(double), "h[i]", comm );
workspace->v[i] = (double*) scalloc( total_cap, sizeof(double), "v[i]", comm );
}
/* CG storage */
workspace->r = (real*) scalloc( total_cap, sizeof(real), "r", comm );
workspace->d = (real*) scalloc( total_cap, sizeof(real), "d", comm );
workspace->q = (real*) scalloc( total_cap, sizeof(real), "q", comm );
workspace->p = (real*) scalloc( total_cap, sizeof(real), "p", comm );
workspace->r = (double*) scalloc( total_cap, sizeof(double), "r", comm );
workspace->d = (double*) scalloc( total_cap, sizeof(double), "d", comm );
workspace->q = (double*) scalloc( total_cap, sizeof(double), "q", comm );
workspace->p = (double*) scalloc( total_cap, sizeof(double), "p", comm );
workspace->r2 = (rvec2*) scalloc( total_cap, sizeof(rvec2), "r2", comm );
workspace->d2 = (rvec2*) scalloc( total_cap, sizeof(rvec2), "d2", comm );
workspace->q2 = (rvec2*) scalloc( total_cap, sizeof(rvec2), "q2", comm );
@ -274,8 +274,8 @@ int Allocate_Workspace( reax_system *system, control_params *control,
// /* force related storage */
workspace->f = (rvec*) scalloc( total_cap, sizeof(rvec), "f", comm );
workspace->CdDelta = (real*)
scalloc( total_cap, sizeof(real), "CdDelta", comm );
workspace->CdDelta = (double*)
scalloc( total_cap, sizeof(double), "CdDelta", comm );
return SUCCESS;
}

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

@ -257,14 +257,14 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
}
int BOp( storage *workspace, reax_list *bonds, real bo_cut,
int BOp( storage *workspace, reax_list *bonds, double bo_cut,
int i, int btop_i, far_neighbor_data *nbr_pj,
single_body_parameters *sbp_i, single_body_parameters *sbp_j,
two_body_parameters *twbp ) {
int j, btop_j;
real r2, C12, C34, C56;
real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
real BO, BO_s, BO_pi, BO_pi2;
double r2, C12, C34, C56;
double Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
double BO, BO_s, BO_pi, BO_pi2;
bond_data *ibond, *jbond;
bond_order_data *bo_ij, *bo_ji;
@ -370,14 +370,14 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
{
int i, j, pj, type_i, type_j;
int start_i, end_i, sym_index;
real val_i, Deltap_i, Deltap_boc_i;
real val_j, Deltap_j, Deltap_boc_j;
real f1, f2, f3, f4, f5, f4f5, exp_f4, exp_f5;
real exp_p1i, exp_p2i, exp_p1j, exp_p2j;
real temp, u1_ij, u1_ji, Cf1A_ij, Cf1B_ij, Cf1_ij, Cf1_ji;
real Cf45_ij, Cf45_ji, p_lp1; //u_ij, u_ji
real A0_ij, A1_ij, A2_ij, A2_ji, A3_ij, A3_ji;
real explp1, p_boc1, p_boc2;
double val_i, Deltap_i, Deltap_boc_i;
double val_j, Deltap_j, Deltap_boc_j;
double f1, f2, f3, f4, f5, f4f5, exp_f4, exp_f5;
double exp_p1i, exp_p2i, exp_p1j, exp_p2j;
double temp, u1_ij, u1_ji, Cf1A_ij, Cf1B_ij, Cf1_ij, Cf1_ji;
double Cf45_ij, Cf45_ji, p_lp1; //u_ij, u_ji
double A0_ij, A1_ij, A2_ij, A2_ji, A3_ij, A3_ji;
double explp1, p_boc1, p_boc2;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
bond_order_data *bo_ij, *bo_ji;

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

@ -30,33 +30,33 @@
#include "reaxc_types.h"
typedef struct{
real C1dbo, C2dbo, C3dbo;
real C1dbopi, C2dbopi, C3dbopi, C4dbopi;
real C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2;
real C1dDelta, C2dDelta, C3dDelta;
double C1dbo, C2dbo, C3dbo;
double C1dbopi, C2dbopi, C3dbopi, C4dbopi;
double C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2;
double C1dDelta, C2dDelta, C3dDelta;
}dbond_coefficients;
#ifdef TEST_FORCES
void Get_dBO( reax_system*, reax_list**, int, int, real, rvec* );
void Get_dBO( reax_system*, reax_list**, int, int, double, rvec* );
void Get_dBOpinpi2( reax_system*, reax_list**,
int, int, real, real, rvec*, rvec* );
int, int, double, double, rvec*, rvec* );
void Add_dBO( reax_system*, reax_list**, int, int, real, rvec* );
void Add_dBO( reax_system*, reax_list**, int, int, double, rvec* );
void Add_dBOpinpi2( reax_system*, reax_list**,
int, int, real, real, rvec*, rvec* );
int, int, double, double, rvec*, rvec* );
void Add_dBO_to_Forces( reax_system*, reax_list**, int, int, real );
void Add_dBO_to_Forces( reax_system*, reax_list**, int, int, double );
void Add_dBOpinpi2_to_Forces( reax_system*, reax_list**,
int, int, real, real );
int, int, double, double );
void Add_dDelta( reax_system*, reax_list**, int, real, rvec* );
void Add_dDelta_to_Forces( reax_system *, reax_list**, int, real );
void Add_dDelta( reax_system*, reax_list**, int, double, rvec* );
void Add_dDelta_to_Forces( reax_system *, reax_list**, int, double );
#endif
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*,
int BOp(storage*, reax_list*, double, int, int, far_neighbor_data*,
single_body_parameters*, single_body_parameters*, two_body_parameters*);
void BO( reax_system*, control_params*, simulation_data*,
storage*, reax_list**, output_controls* );

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

@ -38,10 +38,10 @@ void Bonds( reax_system *system, control_params *control,
int i, j, pj, natoms;
int start_i, end_i;
int type_i, type_j;
real ebond, pow_BOs_be2, exp_be12, CEbo;
real gp3, gp4, gp7, gp10, gp37;
real exphu, exphua1, exphub1, exphuov, hulpov, estriph;
real decobdbo, decobdboua, decobdboub;
double ebond, pow_BOs_be2, exp_be12, CEbo;
double gp3, gp4, gp7, gp10, gp37;
double exphu, exphua1, exphub1, exphuov, hulpov, estriph;
double decobdbo, decobdboua, decobdboub;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
bond_order_data *bo_ij;

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

@ -34,7 +34,7 @@ char Read_Control_File( char *control_file, control_params* control,
FILE *fp;
char *s, **tmp;
int i,ival;
real val;
double val;
/* open control file */
if ( (fp = fopen( control_file, "r" ) ) == NULL ) {

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

@ -38,7 +38,7 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
int c, i, j, k, l, m, n, o, p, cnt;
int lgflag = control->lgflag;
int errorflag = 1;
real val;
double val;
MPI_Comm comm;
comm = MPI_COMM_WORLD;
@ -64,14 +64,14 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
}
reax->gp.n_global = n;
reax->gp.l = (real*) malloc(sizeof(real)*n);
reax->gp.l = (double*) malloc(sizeof(double)*n);
/* see reax_types.h for mapping between l[i] and the lambdas used in ff */
for (i=0; i < n; i++) {
fgets(s,MAX_LINE,fp);
c = Tokenize(s,&tmp);
val = (real) atof(tmp[0]);
val = (double) atof(tmp[0]);
reax->gp.l[i] = val;
}

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

@ -171,9 +171,9 @@ void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists,
}
real Compute_H( real r, real gamma, real *ctap )
double Compute_H( double r, double gamma, double *ctap )
{
real taper, dr3gamij_1, dr3gamij_3;
double taper, dr3gamij_1, dr3gamij_3;
taper = ctap[7] * r + ctap[6];
taper = taper * r + ctap[5];
@ -189,10 +189,10 @@ real Compute_H( real r, real gamma, real *ctap )
}
real Compute_tabH( real r_ij, int ti, int tj )
double Compute_tabH( double r_ij, int ti, int tj )
{
int r, tmin, tmax;
real val, dif, base;
double val, dif, base;
LR_lookup_table *t;
tmin = MIN( ti, tj );
@ -202,7 +202,7 @@ real Compute_tabH( real r_ij, int ti, int tj )
/* cubic spline interpolation */
r = (int)(r_ij * t->inv_dx);
if( r == 0 ) ++r;
base = (real)(r+1) * t->dx;
base = (double)(r+1) * t->dx;
dif = r_ij - base;
val = ((t->ele[r].d*dif + t->ele[r].c)*dif + t->ele[r].b)*dif +
t->ele[r].a;
@ -221,7 +221,7 @@ void Init_Forces( reax_system *system, control_params *control,
int Htop, btop_i, btop_j, num_bonds, num_hbonds;
int ihb, jhb, ihb_top, jhb_top;
int local, flag, renbr;
real r_ij, cutoff;
double r_ij, cutoff;
sparse_matrix *H;
reax_list *far_nbrs, *bonds, *hbonds;
single_body_parameters *sbp_i, *sbp_j;
@ -388,7 +388,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
int btop_i, btop_j, num_bonds, num_hbonds;
int ihb, jhb, ihb_top, jhb_top;
int local, flag, renbr;
real cutoff;
double cutoff;
reax_list *far_nbrs, *bonds, *hbonds;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
@ -531,10 +531,10 @@ void Estimate_Storages( reax_system *system, control_params *control,
int type_i, type_j;
int ihb, jhb;
int local;
real cutoff;
real r_ij;
real C12, C34, C56;
real BO, BO_s, BO_pi, BO_pi2;
double cutoff;
double r_ij;
double C12, C34, C56;
double BO, BO_s, BO_pi, BO_pi2;
reax_list *far_nbrs;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;

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

@ -42,8 +42,8 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
int itr, top;
int num_hb_intrs = 0;
ivec rel_jk;
real r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2;
real e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3;
double r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2;
double e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3;
rvec dcos_theta_di, dcos_theta_dj, dcos_theta_dk;
rvec dvec_jk, force, ext_press;
hbond_parameters *hbp;
@ -55,7 +55,7 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
hbond_data *hbond_list;
// tally variables
real fi_tmp[3], fk_tmp[3], delij[3], delkj[3];
double fi_tmp[3], fk_tmp[3], delij[3], delkj[3];
bonds = (*lists) + BONDS;
bond_list = bonds->select.bond_list;

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

@ -82,9 +82,9 @@ int Init_Simulation_Data( reax_system *system, control_params *control,
void Init_Taper( control_params *control, storage *workspace, MPI_Comm comm )
{
real d1, d7;
real swa, swa2, swa3;
real swb, swb2, swb3;
double d1, d7;
double swa, swa2, swa3;
double swb, swb2, swb3;
swa = control->nonb_low;
swb = control->nonb_cut;

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

@ -430,7 +430,7 @@ void Print_Linear_System( reax_system *system, control_params *control,
}
void Print_LinSys_Soln( reax_system *system, real *x, real *b_prm, real *b )
void Print_LinSys_Soln( reax_system *system, double *x, double *b_prm, double *b )
{
int i;
char fname[100];

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

@ -49,7 +49,7 @@ void Print_Far_Neighbors( reax_system*, reax_list**, control_params *);
void Print_Sparse_Matrix( reax_system*, sparse_matrix* );
void Print_Sparse_Matrix2( reax_system*, sparse_matrix*, char* );
void Print_Linear_System( reax_system*, control_params*, storage*, int );
void Print_LinSys_Soln( reax_system*, real*, real*, real* );
void Print_LinSys_Soln( reax_system*, double*, double*, double* );
void Print_Charges( reax_system* );
void Print_Bonds( reax_system*, reax_list*, char* );
void Print_Bond_List2( reax_system*, reax_list*, char* );

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

@ -31,10 +31,10 @@
LR_lookup_table **LR;
void Tridiagonal_Solve( const real *a, const real *b,
real *c, real *d, real *x, unsigned int n){
void Tridiagonal_Solve( const double *a, const double *b,
double *c, double *d, double *x, unsigned int n){
int i;
real id;
double id;
c[0] /= b[0]; /* Division by zero risk. */
d[0] /= b[0]; /* Division by zero would imply a singular matrix. */
@ -50,19 +50,19 @@ void Tridiagonal_Solve( const real *a, const real *b,
}
void Natural_Cubic_Spline( const real *h, const real *f,
void Natural_Cubic_Spline( const double *h, const double *f,
cubic_spline_coef *coef, unsigned int n,
MPI_Comm comm )
{
int i;
real *a, *b, *c, *d, *v;
double *a, *b, *c, *d, *v;
/* allocate space for the linear system */
a = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm );
b = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm );
c = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm );
d = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm );
v = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm );
a = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
b = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
c = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
d = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
v = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
/* build the linear system */
a[0] = a[1] = a[n-1] = 0;
@ -101,19 +101,19 @@ void Natural_Cubic_Spline( const real *h, const real *f,
void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast,
void Complete_Cubic_Spline( const double *h, const double *f, double v0, double vlast,
cubic_spline_coef *coef, unsigned int n,
MPI_Comm comm )
{
int i;
real *a, *b, *c, *d, *v;
double *a, *b, *c, *d, *v;
/* allocate space for the linear system */
a = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm );
b = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm );
c = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm );
d = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm );
v = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm );
a = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
b = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
c = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
d = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
v = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
/* build the linear system */
a[0] = 0;
@ -150,14 +150,14 @@ void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast,
}
void LR_Lookup( LR_lookup_table *t, real r, LR_data *y )
void LR_Lookup( LR_lookup_table *t, double r, LR_data *y )
{
int i;
real base, dif;
double base, dif;
i = (int)(r * t->inv_dx);
if( i == 0 ) ++i;
base = (real)(i+1) * t->dx;
base = (double)(i+1) * t->dx;
dif = r - base;
y->e_vdW = ((t->vdW[i].d*dif + t->vdW[i].c)*dif + t->vdW[i].b)*dif +
@ -180,9 +180,9 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
int i, j, r;
int num_atom_types;
int existing_types[MAX_ATOM_TYPES], aggregated[MAX_ATOM_TYPES];
real dr;
real *h, *fh, *fvdw, *fele, *fCEvd, *fCEclmb;
real v0_vdw, v0_ele, vlast_vdw, vlast_ele;
double dr;
double *h, *fh, *fvdw, *fele, *fCEvd, *fCEclmb;
double v0_vdw, v0_ele, vlast_vdw, vlast_ele;
MPI_Comm comm;
/* initializations */
@ -194,18 +194,18 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
num_atom_types = system->reax_param.num_atom_types;
dr = control->nonb_cut / control->tabulate;
h = (real*)
smalloc( (control->tabulate+2) * sizeof(real), "lookup:h", comm );
fh = (real*)
smalloc( (control->tabulate+2) * sizeof(real), "lookup:fh", comm );
fvdw = (real*)
smalloc( (control->tabulate+2) * sizeof(real), "lookup:fvdw", comm );
fCEvd = (real*)
smalloc( (control->tabulate+2) * sizeof(real), "lookup:fCEvd", comm );
fele = (real*)
smalloc( (control->tabulate+2) * sizeof(real), "lookup:fele", comm );
fCEclmb = (real*)
smalloc( (control->tabulate+2) * sizeof(real), "lookup:fCEclmb", comm );
h = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:h", comm );
fh = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fh", comm );
fvdw = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fvdw", comm );
fCEvd = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fCEvd", comm );
fele = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fele", comm );
fCEclmb = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fCEclmb", comm );
LR = (LR_lookup_table**)
scalloc( num_atom_types, sizeof(LR_lookup_table*), "lookup:LR", comm );

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

@ -35,17 +35,17 @@ void Atom_Energy( reax_system *system, control_params *control,
output_controls *out_control )
{
int i, j, pj, type_i, type_j;
real Delta_lpcorr, dfvl;
real e_lp, expvd2, inv_expvd2, dElp, CElp, DlpVi;
real e_lph, Di, vov3, deahu2dbo, deahu2dsbo;
real e_ov, CEover1, CEover2, CEover3, CEover4;
real exp_ovun1, exp_ovun2, sum_ovun1, sum_ovun2;
real exp_ovun2n, exp_ovun6, exp_ovun8;
real inv_exp_ovun1, inv_exp_ovun2, inv_exp_ovun2n, inv_exp_ovun8;
real e_un, CEunder1, CEunder2, CEunder3, CEunder4;
real p_lp2, p_lp3;
real p_ovun2, p_ovun3, p_ovun4, p_ovun5, p_ovun6, p_ovun7, p_ovun8;
real eng_tmp;
double Delta_lpcorr, dfvl;
double e_lp, expvd2, inv_expvd2, dElp, CElp, DlpVi;
double e_lph, Di, vov3, deahu2dbo, deahu2dsbo;
double e_ov, CEover1, CEover2, CEover3, CEover4;
double exp_ovun1, exp_ovun2, sum_ovun1, sum_ovun2;
double exp_ovun2n, exp_ovun6, exp_ovun8;
double inv_exp_ovun1, inv_exp_ovun2, inv_exp_ovun2n, inv_exp_ovun8;
double e_un, CEunder1, CEunder2, CEunder3, CEunder4;
double p_lp2, p_lp3;
double p_ovun2, p_ovun3, p_ovun4, p_ovun5, p_ovun6, p_ovun7, p_ovun8;
double eng_tmp;
int numbonds;
single_body_parameters *sbp_i;

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

@ -38,20 +38,20 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
int i, j, pj, natoms;
int start_i, end_i, flag;
rc_tagint orig_i, orig_j;
real p_vdW1, p_vdW1i;
real powr_vdW1, powgi_vdW1;
real tmp, r_ij, fn13, exp1, exp2;
real Tap, dTap, dfn13, CEvd, CEclmb, de_core;
real dr3gamij_1, dr3gamij_3;
real e_ele, e_vdW, e_core, SMALL = 0.0001;
real e_lg, de_lg, r_ij5, r_ij6, re6;
double p_vdW1, p_vdW1i;
double powr_vdW1, powgi_vdW1;
double tmp, r_ij, fn13, exp1, exp2;
double Tap, dTap, dfn13, CEvd, CEclmb, de_core;
double dr3gamij_1, dr3gamij_3;
double e_ele, e_vdW, e_core, SMALL = 0.0001;
double e_lg, de_lg, r_ij5, r_ij6, re6;
rvec temp, ext_press;
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
reax_list *far_nbrs;
// Tallying variables:
real pe_vdw, f_tmp, delij[3];
double pe_vdw, f_tmp, delij[3];
natoms = system->n;
far_nbrs = (*lists) + FAR_NBRS;
@ -212,10 +212,10 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
int type_i, type_j, tmin, tmax;
int start_i, end_i, flag;
rc_tagint orig_i, orig_j;
real r_ij, base, dif;
real e_vdW, e_ele;
real CEvd, CEclmb, SMALL = 0.0001;
real f_tmp, delij[3];
double r_ij, base, dif;
double e_vdW, e_ele;
double CEvd, CEclmb, SMALL = 0.0001;
double f_tmp, delij[3];
rvec temp, ext_press;
far_neighbor_data *nbr_pj;
@ -265,7 +265,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;
base = (real)(r+1) * t->dx;
base = (double)(r+1) * t->dx;
dif = r_ij - base;
e_vdW = ((t->vdW[r].d*dif + t->vdW[r].c)*dif + t->vdW[r].b)*dif +
@ -319,7 +319,7 @@ 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, en_tmp;
double q, en_tmp;
data->my_en.e_pol = 0.0;
for( i = 0; i < system->n; i++ ) {
@ -338,16 +338,16 @@ void Compute_Polarization_Energy( reax_system *system, simulation_data *data )
}
void LR_vdW_Coulomb( reax_system *system, storage *workspace,
control_params *control, int i, int j, real r_ij, LR_data *lr )
control_params *control, int i, int j, double r_ij, LR_data *lr )
{
real p_vdW1 = system->reax_param.gp.l[28];
real p_vdW1i = 1.0 / p_vdW1;
real powr_vdW1, powgi_vdW1;
real tmp, fn13, exp1, exp2;
real Tap, dTap, dfn13;
real dr3gamij_1, dr3gamij_3;
real e_core, de_core;
real e_lg, de_lg, r_ij5, r_ij6, re6;
double p_vdW1 = system->reax_param.gp.l[28];
double p_vdW1i = 1.0 / p_vdW1;
double powr_vdW1, powgi_vdW1;
double tmp, fn13, exp1, exp2;
double Tap, dTap, dfn13;
double dr3gamij_1, dr3gamij_3;
double e_core, de_core;
double e_lg, de_lg, r_ij5, r_ij6, re6;
two_body_parameters *twbp;
twbp = &(system->reax_param.tbp[i][j]);

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

@ -39,5 +39,5 @@ void Tabulated_vdW_Coulomb_Energy( reax_system*, control_params*,
void Compute_Polarization_Energy( reax_system*, simulation_data* );
void LR_vdW_Coulomb( reax_system*, storage*, control_params*,
int, int, real, LR_data* );
int, int, double, LR_data* );
#endif

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

@ -111,9 +111,9 @@ void Reset_Timing( reax_timing *rt )
void Reset_Workspace( reax_system *system, storage *workspace )
{
memset( workspace->total_bond_order, 0, system->total_cap * sizeof( real ) );
memset( workspace->total_bond_order, 0, system->total_cap * sizeof( double ) );
memset( workspace->dDeltap_self, 0, system->total_cap * sizeof( rvec ) );
memset( workspace->CdDelta, 0, system->total_cap * sizeof( real ) );
memset( workspace->CdDelta, 0, system->total_cap * sizeof( double ) );
memset( workspace->f, 0, system->total_cap * sizeof( rvec ) );
}

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

@ -31,7 +31,7 @@
void Temperature_Control( control_params *control, simulation_data *data )
{
real tmp;
double tmp;
if( control->T_mode == 1 ) {// step-wise temperature control
if((data->step-data->prev_steps) % ((int)(control->T_freq/control->dt))==0){
@ -55,7 +55,7 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data,
{
int i;
rvec p;
real m;
double m;
data->my_en.e_kin = 0.0;
data->sys_en.e_kin = 0.0;
@ -82,7 +82,7 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data,
void Compute_System_Energy( reax_system *system, simulation_data *data,
MPI_Comm comm )
{
real my_en[15], sys_en[15];
double my_en[15], sys_en[15];
my_en[0] = data->my_en.e_bond;
my_en[1] = data->my_en.e_ov;
@ -141,7 +141,7 @@ void Compute_Total_Mass( reax_system *system, simulation_data *data,
MPI_Comm comm )
{
int i;
real tmp;
double tmp;
tmp = 0;
for( i = 0; i < system->n; i++ )
@ -158,8 +158,8 @@ void Compute_Center_of_Mass( reax_system *system, simulation_data *data,
mpi_datatypes *mpi_data, MPI_Comm comm )
{
int i;
real m, det; //xx, xy, xz, yy, yz, zz;
real tmp_mat[6], tot_mat[6];
double m, det; //xx, xy, xz, yy, yz, zz;
double tmp_mat[6], tot_mat[6];
rvec my_xcm, my_vcm, my_amcm, my_avcm;
rvec tvec, diff;
rtensor mat, inv;

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

@ -30,7 +30,7 @@
void Transform( rvec x1, simulation_box *box, char flag, rvec x2 )
{
int i, j;
real tmp;
double tmp;
if (flag > 0) {
for (i=0; i < 3; i++) {
@ -79,16 +79,16 @@ void Fit_to_Periodic_Box( simulation_box *box, rvec *p )
}
struct timeval tim;
real t_end;
double t_end;
real Get_Time( )
double Get_Time( )
{
gettimeofday(&tim, NULL );
return( tim.tv_sec + (tim.tv_usec / 1000000.0) );
}
real Get_Timing_Info( real t_start )
double Get_Timing_Info( double t_start )
{
gettimeofday(&tim, NULL );
t_end = tim.tv_sec + (tim.tv_usec / 1000000.0);
@ -96,7 +96,7 @@ real Get_Timing_Info( real t_start )
}
void Update_Timing_Info( real *t_start, real *timing )
void Update_Timing_Info( double *t_start, double *timing )
{
gettimeofday(&tim, NULL );
t_end = tim.tv_sec + (tim.tv_usec / 1000000.0);

10
src/USER-REAXC/reaxc_tool_box.h
View File

@ -44,14 +44,14 @@ int iown_midpoint( simulation_box*, rvec, rvec );
/* from grid.h */
void GridCell_Closest_Point( grid_cell*, grid_cell*, ivec, ivec, rvec );
void GridCell_to_Box_Points( grid_cell*, ivec, rvec, rvec );
real DistSqr_between_Special_Points( rvec, rvec );
real DistSqr_to_Special_Point( rvec, rvec );
double DistSqr_between_Special_Points( rvec, rvec );
double DistSqr_to_Special_Point( rvec, rvec );
int Relative_Coord_Encoding( ivec );
/* from system_props.h */
real Get_Time( );
real Get_Timing_Info( real );
void Update_Timing_Info( real*, real* );
double Get_Time( );
double Get_Timing_Info( double );
void Update_Timing_Info( double*, double* );
/* from io_tools.h */
int Get_Atom_Type( reax_interaction*, char*, MPI_Comm );

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

@ -33,20 +33,20 @@
#define MIN_SINE 1e-10
real Calculate_Omega( rvec dvec_ij, real r_ij,
rvec dvec_jk, real r_jk,
rvec dvec_kl, real r_kl,
rvec dvec_li, real r_li,
double Calculate_Omega( rvec dvec_ij, double r_ij,
rvec dvec_jk, double r_jk,
rvec dvec_kl, double r_kl,
rvec dvec_li, double r_li,
three_body_interaction_data *p_ijk,
three_body_interaction_data *p_jkl,
rvec dcos_omega_di, rvec dcos_omega_dj,
rvec dcos_omega_dk, rvec dcos_omega_dl,
output_controls *out_control )
{
real unnorm_cos_omega, unnorm_sin_omega, omega;
real sin_ijk, cos_ijk, sin_jkl, cos_jkl;
real htra, htrb, htrc, hthd, hthe, hnra, hnrc, hnhd, hnhe;
real arg, poem, tel;
double unnorm_cos_omega, unnorm_sin_omega, omega;
double sin_ijk, cos_ijk, sin_jkl, cos_jkl;
double htra, htrb, htrc, hthd, hthe, hnra, hnrc, hnhd, hnhe;
double arg, poem, tel;
rvec cross_jk_kl;
sin_ijk = sin( p_ijk->theta );
@ -128,25 +128,25 @@ void Torsion_Angles( reax_system *system, control_params *control,
int start_pj, end_pj, start_pk, end_pk;
int num_frb_intrs = 0;
real Delta_j, Delta_k;
real r_ij, r_jk, r_kl, r_li;
real BOA_ij, BOA_jk, BOA_kl;
double Delta_j, Delta_k;
double r_ij, r_jk, r_kl, r_li;
double BOA_ij, BOA_jk, BOA_kl;
real exp_tor2_ij, exp_tor2_jk, exp_tor2_kl;
real exp_tor1, exp_tor3_DjDk, exp_tor4_DjDk, exp_tor34_inv;
real exp_cot2_jk, exp_cot2_ij, exp_cot2_kl;
real fn10, f11_DjDk, dfn11, fn12;
real theta_ijk, theta_jkl;
real sin_ijk, sin_jkl;
real cos_ijk, cos_jkl;
real tan_ijk_i, tan_jkl_i;
real omega, cos_omega, cos2omega, cos3omega;
double exp_tor2_ij, exp_tor2_jk, exp_tor2_kl;
double exp_tor1, exp_tor3_DjDk, exp_tor4_DjDk, exp_tor34_inv;
double exp_cot2_jk, exp_cot2_ij, exp_cot2_kl;
double fn10, f11_DjDk, dfn11, fn12;
double theta_ijk, theta_jkl;
double sin_ijk, sin_jkl;
double cos_ijk, cos_jkl;
double tan_ijk_i, tan_jkl_i;
double omega, cos_omega, cos2omega, cos3omega;
rvec dcos_omega_di, dcos_omega_dj, dcos_omega_dk, dcos_omega_dl;
real CV, cmn, CEtors1, CEtors2, CEtors3, CEtors4;
real CEtors5, CEtors6, CEtors7, CEtors8, CEtors9;
real Cconj, CEconj1, CEconj2, CEconj3;
real CEconj4, CEconj5, CEconj6;
real e_tor, e_con;
double CV, cmn, CEtors1, CEtors2, CEtors3, CEtors4;
double CEtors5, CEtors6, CEtors7, CEtors8, CEtors9;
double Cconj, CEconj1, CEconj2, CEconj3;
double CEconj4, CEconj5, CEconj6;
double e_tor, e_con;
rvec dvec_li;
rvec force, ext_press;
ivec rel_box_jl;
@ -156,16 +156,16 @@ void Torsion_Angles( reax_system *system, control_params *control,
bond_data *pbond_ij, *pbond_jk, *pbond_kl;
bond_order_data *bo_ij, *bo_jk, *bo_kl;
three_body_interaction_data *p_ijk, *p_jkl;
real p_tor2 = system->reax_param.gp.l[23];
real p_tor3 = system->reax_param.gp.l[24];
real p_tor4 = system->reax_param.gp.l[25];
real p_cot2 = system->reax_param.gp.l[27];
double p_tor2 = system->reax_param.gp.l[23];
double p_tor3 = system->reax_param.gp.l[24];
double p_tor4 = system->reax_param.gp.l[25];
double p_cot2 = system->reax_param.gp.l[27];
reax_list *bonds = (*lists) + BONDS;
reax_list *thb_intrs = (*lists) + THREE_BODIES;
// Virial tallying variables
real delil[3], deljl[3], delkl[3];
real eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3];
double delil[3], deljl[3], delkl[3];
double eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3];
natoms = system->n;

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

@ -62,11 +62,10 @@
/********************** TYPE DEFINITIONS ********************/
typedef int ivec[3];
typedef double real;
typedef real rvec[3];
typedef real rtensor[3][3];
typedef real rvec2[2];
typedef real rvec4[4];
typedef double rvec[3];
typedef double rtensor[3][3];
typedef double rvec2[2];
typedef double rvec4[4];
// import LAMMPS' definition of tagint
@ -111,7 +110,7 @@ typedef struct
typedef struct
{
int n_global;
real* l;
double* l;
int vdw_type;
} global_parameters;
@ -122,44 +121,44 @@ typedef struct
/* Line one in field file */
char name[15]; // Two character atom name
real r_s;
real valency; // Valency of the atom
real mass; // Mass of atom
real r_vdw;
real epsilon;
real gamma;
real r_pi;
real valency_e;
real nlp_opt;
double r_s;
double valency; // Valency of the atom
double mass; // Mass of atom
double r_vdw;
double epsilon;
double gamma;
double r_pi;
double valency_e;
double nlp_opt;
/* Line two in field file */
real alpha;
real gamma_w;
real valency_boc;
real p_ovun5;
real chi;
real eta;
double alpha;
double gamma_w;
double valency_boc;
double p_ovun5;
double chi;
double eta;
int p_hbond; // 1 for H, 2 for hbonding atoms (O,S,P,N), 0 for others
/* Line three in field file */
real r_pi_pi;
real p_lp2;
real b_o_131;
real b_o_132;
real b_o_133;
double r_pi_pi;
double p_lp2;
double b_o_131;
double b_o_132;
double b_o_133;
/* Line four in the field file */
real p_ovun2;
real p_val3;
real valency_val;
real p_val5;
real rcore2;
real ecore2;
real acore2;
double p_ovun2;
double p_val3;
double valency_val;
double p_val5;
double rcore2;
double ecore2;
double acore2;
/* Line five in the ffield file, only for lgvdw yes */
real lgcij;
real lgre;
double lgcij;
double lgre;
} single_body_parameters;
@ -168,29 +167,29 @@ typedef struct
/* Two Body Parameters */
typedef struct {
/* Bond Order parameters */
real p_bo1,p_bo2,p_bo3,p_bo4,p_bo5,p_bo6;
real r_s, r_p, r_pp; // r_o distances in BO formula
real p_boc3, p_boc4, p_boc5;
double p_bo1,p_bo2,p_bo3,p_bo4,p_bo5,p_bo6;
double r_s, r_p, r_pp; // r_o distances in BO formula
double p_boc3, p_boc4, p_boc5;
/* Bond Energy parameters */
real p_be1, p_be2;
real De_s, De_p, De_pp;
double p_be1, p_be2;
double De_s, De_p, De_pp;
/* Over/Under coordination parameters */
real p_ovun1;
double p_ovun1;
/* Van der Waal interaction parameters */
real D;
real alpha;
real r_vdW;
real gamma_w;
real rcore, ecore, acore;
real lgcij, lgre;
double D;
double alpha;
double r_vdW;
double gamma_w;
double rcore, ecore, acore;
double lgcij, lgre;
/* electrostatic parameters */
real gamma; // note: this parameter is gamma^-3 and not gamma.
double gamma; // note: this parameter is gamma^-3 and not gamma.
real v13cor, ovc;
double v13cor, ovc;
} two_body_parameters;
@ -198,14 +197,14 @@ typedef struct {
/* 3-body parameters */
typedef struct {
/* valence angle */
real theta_00;
real p_val1, p_val2, p_val4, p_val7;
double theta_00;
double p_val1, p_val2, p_val4, p_val7;
/* penalty */
real p_pen1;
double p_pen1;
/* 3-body conjugation */
real p_coa1;
double p_coa1;
} three_body_parameters;
@ -218,20 +217,20 @@ typedef struct{
/* hydrogen-bond parameters */
typedef struct{
real r0_hb, p_hb1, p_hb2, p_hb3;
double r0_hb, p_hb1, p_hb2, p_hb3;
} hbond_parameters;
/* 4-body parameters */
typedef struct {
real V1, V2, V3;
double V1, V2, V3;
/* torsion angle */
real p_tor1;
double p_tor1;
/* 4-body conjugation */
real p_cot1;
double p_cot1;
} four_body_parameters;
@ -267,7 +266,7 @@ struct _reax_atom
rvec f; // force
rvec f_old;
real q; // charge
double q; // charge
rvec4 s; // they take part in
rvec4 t; // computing q
@ -287,7 +286,7 @@ typedef _reax_atom reax_atom;
typedef struct
{
real V;
double V;
rvec min, max, box_norms;
rtensor box, box_inv;
@ -299,7 +298,7 @@ typedef struct
struct grid_cell
{
real cutoff;
double cutoff;
rvec min, max;
ivec rel_box;
@ -332,7 +331,7 @@ typedef struct
ivec native_str;
ivec native_end;
real ghost_cut;
double ghost_cut;
ivec ghost_span;
ivec ghost_nonb_span;
ivec ghost_hbond_span;
@ -372,10 +371,10 @@ typedef struct
typedef struct
{
real ghost_nonb;
real ghost_hbond;
real ghost_bond;
real ghost_cutoff;
double ghost_nonb;
double ghost_hbond;
double ghost_bond;
double ghost_cutoff;
} boundary_cutoff;
using LAMMPS_NS::Pair;
@ -399,7 +398,7 @@ struct _reax_system
class Pair *pair_ptr;
int my_bonds;
int mincap;
real safezone, saferzone;
double safezone, saferzone;
};
typedef _reax_system reax_system;
@ -421,7 +420,7 @@ typedef struct
5 : NPT (Parrinello-Rehman-Nose-Hoover) Anisotropic*/
int ensemble;
int nsteps;
real dt;
double dt;
int geo_format;
int restart;
@ -431,33 +430,33 @@ typedef struct
int reposition_atoms;
int reneighbor;
real vlist_cut;
real bond_cut;
real nonb_cut, nonb_low;
real hbond_cut;
real user_ghost_cut;
double vlist_cut;
double bond_cut;
double nonb_cut, nonb_low;
double hbond_cut;
double user_ghost_cut;
real bg_cut;
real bo_cut;
real thb_cut;
real thb_cutsq;
double bg_cut;
double bo_cut;
double thb_cut;
double thb_cutsq;
int tabulate;
int qeq_freq;
real q_err;
double q_err;
int refactor;
real droptol;
double droptol;
real T_init, T_final, T;
real Tau_T;
double T_init, T_final, T;
double Tau_T;
int T_mode;
real T_rate, T_freq;
double T_rate, T_freq;
int virial;
rvec P, Tau_P, Tau_PT;
int press_mode;
real compressibility;
double compressibility;
int molecular_analysis;
int num_ignored;
@ -476,22 +475,22 @@ typedef struct
typedef struct
{
real T;
real xi;
real v_xi;
real v_xi_old;
real G_xi;
double T;
double xi;
double v_xi;
double v_xi_old;
double G_xi;
} thermostat;
typedef struct
{
real P;
real eps;
real v_eps;
real v_eps_old;
real a_eps;
double P;
double eps;
double v_eps;
double v_eps_old;
double a_eps;
} isotropic_barostat;
@ -499,12 +498,12 @@ typedef struct
typedef struct
{
rtensor P;
real P_scalar;
double P_scalar;
real eps;
real v_eps;
real v_eps_old;
real a_eps;
double eps;
double v_eps;
double v_eps_old;
double a_eps;
rtensor h0;
rtensor v_g0;
@ -516,17 +515,17 @@ typedef struct
typedef struct
{
real start;
real end;
real elapsed;
double start;
double end;
double elapsed;
real total;
real comm;
real nbrs;
real init_forces;
real bonded;
real nonb;
real qEq;
double total;
double comm;
double nbrs;
double init_forces;
double bonded;
double nonb;
double qEq;
int s_matvecs;
int t_matvecs;
} reax_timing;
@ -534,41 +533,41 @@ typedef struct
typedef struct
{
real e_tot;
real e_kin; // Total kinetic energy
real e_pot;
double e_tot;
double e_kin; // Total kinetic energy
double e_pot;
real e_bond; // Total bond energy
real e_ov; // Total over coordination
real e_un; // Total under coordination energy
real e_lp; // Total under coordination energy
real e_ang; // Total valance angle energy
real e_pen; // Total penalty energy
real e_coa; // Total three body conjgation energy
real e_hb; // Total Hydrogen bond energy
real e_tor; // Total torsional energy
real e_con; // Total four body conjugation energy
real e_vdW; // Total van der Waals energy
real e_ele; // Total electrostatics energy
real e_pol; // Polarization energy
double e_bond; // Total bond energy
double e_ov; // Total over coordination
double e_un; // Total under coordination energy
double e_lp; // Total under coordination energy
double e_ang; // Total valance angle energy
double e_pen; // Total penalty energy
double e_coa; // Total three body conjgation energy
double e_hb; // Total Hydrogen bond energy
double e_tor; // Total torsional energy
double e_con; // Total four body conjugation energy
double e_vdW; // Total van der Waals energy
double e_ele; // Total electrostatics energy
double e_pol; // Polarization energy
} energy_data;
typedef struct
{
int step;
int prev_steps;
real time;
double time;
real M; // Total Mass
real inv_M; // 1 / Total Mass
double M; // Total Mass
double inv_M; // 1 / Total Mass
rvec xcm; // Center of mass
rvec vcm; // Center of mass velocity
rvec fcm; // Center of mass force
rvec amcm; // Angular momentum of CoM
rvec avcm; // Angular velocity of CoM
real etran_cm; // Translational kinetic energy of CoM
real erot_cm; // Rotational kinetic energy of CoM
double etran_cm; // Translational kinetic energy of CoM
double erot_cm; // Rotational kinetic energy of CoM
rtensor kinetic; // Kinetic energy tensor
rtensor virial; // Hydrodynamic virial
@ -576,15 +575,15 @@ typedef struct
energy_data my_en;
energy_data sys_en;
real N_f; //Number of degrees of freedom
double N_f; //Number of degrees of freedom
rvec t_scale;
rtensor p_scale;
thermostat therm; // Used in Nose_Hoover method
isotropic_barostat iso_bar;
flexible_barostat flex_bar;
real inv_W;
double inv_W;
real kin_press;
double kin_press;
rvec int_press;
rvec my_ext_press;
rvec ext_press;
@ -597,7 +596,7 @@ typedef struct
typedef struct{
int thb;
int pthb; // pointer to the third body on the central atom's nbrlist
real theta, cos_theta;
double theta, cos_theta;
rvec dcos_di, dcos_dj, dcos_dk;
} three_body_interaction_data;
@ -605,7 +604,7 @@ typedef struct{
typedef struct {
int nbr;
ivec rel_box;
real d;
double d;
rvec dvec;
} far_neighbor_data;
@ -629,11 +628,11 @@ typedef struct{
} dbond_data;
typedef struct{
real BO, BO_s, BO_pi, BO_pi2;
real Cdbo, Cdbopi, Cdbopi2;
real C1dbo, C2dbo, C3dbo;
real C1dbopi, C2dbopi, C3dbopi, C4dbopi;
real C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2;
double BO, BO_s, BO_pi, BO_pi2;
double Cdbo, Cdbopi, Cdbopi2;
double C1dbo, C2dbo, C3dbo;
double C1dbopi, C2dbopi, C3dbopi, C4dbopi;
double C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2;
rvec dBOp, dln_BOp_s, dln_BOp_pi, dln_BOp_pi2;
} bond_order_data;
@ -643,7 +642,7 @@ typedef struct {
int dbond_index;
ivec rel_box;
// rvec ext_factor;
real d;
double d;
rvec dvec;
bond_order_data bo_data;
} bond_data;
@ -651,7 +650,7 @@ typedef struct {
typedef struct {
int j;
real val;
double val;
} sparse_matrix_entry;
typedef struct {
@ -676,35 +675,35 @@ typedef struct
int allocated;
/* communication storage */
real *tmp_dbl[REAX_MAX_NBRS];
double *tmp_dbl[REAX_MAX_NBRS];
rvec *tmp_rvec[REAX_MAX_NBRS];
rvec2 *tmp_rvec2[REAX_MAX_NBRS];
int *within_bond_box;
/* bond order related storage */
real *total_bond_order;
real *Deltap, *Deltap_boc;
real *Delta, *Delta_lp, *Delta_lp_temp, *Delta_e, *Delta_boc, *Delta_val;
real *dDelta_lp, *dDelta_lp_temp;
real *nlp, *nlp_temp, *Clp, *vlpex;
double *total_bond_order;
double *Deltap, *Deltap_boc;
double *Delta, *Delta_lp, *Delta_lp_temp, *Delta_e, *Delta_boc, *Delta_val;
double *dDelta_lp, *dDelta_lp_temp;
double *nlp, *nlp_temp, *Clp, *vlpex;
rvec *dDeltap_self;
int *bond_mark, *done_after;
/* QEq storage */
sparse_matrix *H, *L, *U;
real *Hdia_inv, *b_s, *b_t, *b_prc, *b_prm, *s, *t;
real *droptol;
double *Hdia_inv, *b_s, *b_t, *b_prc, *b_prm, *s, *t;
double *droptol;
rvec2 *b, *x;
/* GMRES storage */
real *y, *z, *g;
real *hc, *hs;
real **h, **v;
double *y, *z, *g;
double *hc, *hs;
double **h, **v;
/* CG storage */
real *r, *d, *q, *p;
double *r, *d, *q, *p;
rvec2 *r2, *d2, *q2, *p2;
/* Taper */
real Tap[8]; //Tap7, Tap6, Tap5, Tap4, Tap3, Tap2, Tap1, Tap0;
double Tap[8]; //Tap7, Tap6, Tap5, Tap4, Tap3, Tap2, Tap1, Tap0;
/* storage for analysis */
int *mark, *old_mark;
@ -718,7 +717,7 @@ typedef struct
rvec *v_const;
/* force calculations */
real *CdDelta; // coefficient of dDelta
double *CdDelta; // coefficient of dDelta
rvec *f;
reallocate_data realloc;
@ -802,27 +801,27 @@ typedef struct
typedef struct
{
real H;
real e_vdW, CEvd;
real e_ele, CEclmb;
double H;
double e_vdW, CEvd;
double e_ele, CEclmb;
} LR_data;
typedef struct
{
real a, b, c, d;
double a, b, c, d;
} cubic_spline_coef;
typedef struct
{
real xmin, xmax;
double xmin, xmax;
int n;
real dx, inv_dx;
real a;
real m;
real c;
double dx, inv_dx;
double a;
double m;
double c;
LR_data *y;
cubic_spline_coef *H;
@ -844,7 +843,7 @@ typedef void (*print_interaction)(reax_system*, control_params*,
simulation_data*, storage*,
reax_list**, output_controls*);
typedef real (*lookup_function)(real);
typedef double (*lookup_function)(double);
typedef void (*message_sorter) (reax_system*, int, int, int, mpi_out_data*);
typedef void (*unpacker) ( reax_system*, int, void*, int, neighbor_proc*, int );

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

@ -30,9 +30,9 @@
#include "reaxc_list.h"
#include "reaxc_vector.h"
static real Dot( real* v1, real* v2, int k )
static double Dot( double* v1, double* v2, int k )
{
real ret = 0.0;
double ret = 0.0;
for( int i=0; i < k; ++i )
ret += v1[i] * v2[i];
@ -40,8 +40,8 @@ static real Dot( real* v1, real* v2, int k )
return ret;
}
void Calculate_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
real *theta, real *cos_theta )
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;
@ -50,18 +50,18 @@ void Calculate_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
(*theta) = acos( *cos_theta );
}
void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
void Calculate_dCos_Theta( rvec dvec_ji, double d_ji, rvec dvec_jk, double 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.0 );
real dot_dvecs = Dot( dvec_ji, dvec_jk, 3 );
real Cdot_inv3 = dot_dvecs * inv_dists3;
double sqr_d_ji = SQR(d_ji);
double sqr_d_jk = SQR(d_jk);
double inv_dists = 1.0 / (d_ji * d_jk);
double inv_dists3 = pow( inv_dists, 3.0 );
double dot_dvecs = Dot( dvec_ji, dvec_jk, 3 );
double Cdot_inv3 = dot_dvecs * inv_dists3;
for( t = 0; t < 3; ++t ) {
(*dcos_theta_di)[t] = dvec_jk[t] * inv_dists -
@ -83,27 +83,27 @@ void Valence_Angles( reax_system *system, control_params *control,
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 BOA_ij, BOA_jk;
double temp, temp_bo_jt, pBOjt7;
double p_val1, p_val2, p_val3, p_val4, p_val5;
double p_val6, p_val7, p_val8, p_val9, p_val10;
double p_pen1, p_pen2, p_pen3, p_pen4;
double p_coa1, p_coa2, p_coa3, p_coa4;
double trm8, expval6, expval7, expval2theta, expval12theta, exp3ij, exp3jk;
double exp_pen2ij, exp_pen2jk, exp_pen3, exp_pen4, trm_pen34, exp_coa2;
double dSBO1, dSBO2, SBO, SBO2, CSBO2, SBOp, prod_SBO, vlpadj;
double CEval1, CEval2, CEval3, CEval4, CEval5, CEval6, CEval7, CEval8;
double CEpen1, CEpen2, CEpen3;
double e_ang, e_coa, e_pen;
double CEcoa1, CEcoa2, CEcoa3, CEcoa4, CEcoa5;
double Cf7ij, Cf7jk, Cf8j, Cf9j;
double f7_ij, f7_jk, f8_Dj, f9_Dj;
double Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta;
double BOA_ij, BOA_jk;
rvec force, ext_press;
// Tallying variables
real eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3];
real delij[3], delkj[3];
double eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3];
double delij[3], delkj[3];
three_body_header *thbh;
three_body_parameters *thbp;

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

@ -32,8 +32,8 @@
void Valence_Angles( reax_system*, control_params*, simulation_data*,
storage*, reax_list**, output_controls* );
void Calculate_Theta( rvec, real, rvec, real, real*, real* );
void Calculate_Theta( rvec, double, rvec, double, double*, double* );
void Calculate_dCos_Theta( rvec, real, rvec, real, rvec*, rvec*, rvec* );
void Calculate_dCos_Theta( rvec, double, rvec, double, rvec*, rvec*, rvec* );
#endif

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

@ -34,7 +34,7 @@ void rvec_Copy( rvec dest, rvec src )
}
void rvec_Scale( rvec ret, real c, rvec v )
void rvec_Scale( rvec ret, double c, rvec v )
{
ret[0] = c * v[0], ret[1] = c * v[1], ret[2] = c * v[2];
}
@ -46,7 +46,7 @@ void rvec_Add( rvec ret, rvec v )
}
void rvec_ScaledAdd( rvec ret, real c, rvec v )
void rvec_ScaledAdd( rvec ret, double c, rvec v )
{
ret[0] += c * v[0], ret[1] += c * v[1], ret[2] += c * v[2];
}
@ -60,7 +60,7 @@ void rvec_Sum( rvec ret, rvec v1 ,rvec v2 )
}
void rvec_ScaledSum( rvec ret, real c1, rvec v1 ,real c2, rvec v2 )
void rvec_ScaledSum( rvec ret, double c1, rvec v1 ,double c2, rvec v2 )
{
ret[0] = c1 * v1[0] + c2 * v2[0];
ret[1] = c1 * v1[1] + c2 * v2[1];
@ -68,13 +68,13 @@ void rvec_ScaledSum( rvec ret, real c1, rvec v1 ,real c2, rvec v2 )
}
real rvec_Dot( rvec v1, rvec v2 )
double rvec_Dot( rvec v1, rvec v2 )
{
return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}
real rvec_ScaledDot( real c1, rvec v1, real c2, rvec v2 )
double rvec_ScaledDot( double c1, rvec v1, double c2, rvec v2 )
{
return (c1*c2) * (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]);
}
@ -138,13 +138,13 @@ void rvec_OuterProduct( rtensor r, rvec v1, rvec v2 )
}
real rvec_Norm_Sqr( rvec v )
double rvec_Norm_Sqr( rvec v )
{
return SQR(v[0]) + SQR(v[1]) + SQR(v[2]);
}
real rvec_Norm( rvec v )
double rvec_Norm( rvec v )
{
return sqrt( SQR(v[0]) + SQR(v[1]) + SQR(v[2]) );
}
@ -187,7 +187,7 @@ void rtensor_MatVec( rvec ret, rtensor m, rvec v )
}
void rtensor_Scale( rtensor ret, real c, rtensor m )
void rtensor_Scale( rtensor ret, double c, rtensor m )
{
int i, j;
@ -217,7 +217,7 @@ void ivec_Copy( ivec dest, ivec src )
}
void ivec_Scale( ivec dest, real C, ivec src )
void ivec_Scale( ivec dest, double C, ivec src )
{
dest[0] = (int)(C * src[0]);
dest[1] = (int)(C * src[1]);

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

@ -32,13 +32,13 @@
#include "reaxc_defs.h"
void rvec_Copy( rvec, rvec );
void rvec_Scale( rvec, real, rvec );
void rvec_Scale( rvec, double, rvec );
void rvec_Add( rvec, rvec );
void rvec_ScaledAdd( rvec, real, rvec );
void rvec_ScaledAdd( rvec, double, rvec );
void rvec_Sum( rvec, rvec, rvec );
void rvec_ScaledSum( rvec, real, rvec, real, rvec );
real rvec_Dot( rvec, rvec );
real rvec_ScaledDot( real, rvec, real, rvec );
void rvec_ScaledSum( rvec, double, rvec, double, rvec );
double rvec_Dot( rvec, rvec );
double rvec_ScaledDot( double, rvec, double, rvec );
void rvec_Multiply( rvec, rvec, rvec );
void rvec_iMultiply( rvec, ivec, rvec );
void rvec_Divide( rvec, rvec, rvec );
@ -46,19 +46,19 @@ void rvec_iDivide( rvec, rvec, ivec );
void rvec_Invert( rvec, rvec );
void rvec_Cross( rvec, rvec, rvec );
void rvec_OuterProduct( rtensor, rvec, rvec );
real rvec_Norm_Sqr( rvec );
real rvec_Norm( rvec );
double rvec_Norm_Sqr( rvec );
double rvec_Norm( rvec );
int rvec_isZero( rvec );
void rvec_MakeZero( rvec );
void rvec_Random( rvec );
void rtensor_MakeZero( rtensor );
void rtensor_MatVec( rvec, rtensor, rvec );
void rtensor_Scale( rtensor, real, rtensor );
void rtensor_Scale( rtensor, double, rtensor );
void ivec_MakeZero( ivec );
void ivec_Copy( ivec, ivec );
void ivec_Scale( ivec, real, ivec );
void ivec_Scale( ivec, double, ivec );
void ivec_Sum( ivec, ivec, ivec );
#endif

14
src/USER-SMD/compute_smd_tlsph_strain_rate.cpp
View File

@ -37,6 +37,7 @@
using namespace Eigen;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHStrainRate::ComputeSMDTLSPHStrainRate(LAMMPS *lmp, int narg, char **arg) :
@ -73,7 +74,6 @@ void ComputeSMDTLSPHStrainRate::init() {
void ComputeSMDTLSPHStrainRate::compute_peratom() {
invoked_peratom = update->ntimestep;
int *mol = atom->molecule;
// grow vector array if necessary
@ -84,8 +84,6 @@ void ComputeSMDTLSPHStrainRate::compute_peratom() {
array_atom = strain_rate_array;
}
int ulsph_flag = 0;
int tlsph_flag = 0;
int itmp = 0;
Matrix3d *D = (Matrix3d *) force->pair->extract("smd/tlsph/strain_rate_ptr", itmp);
if (D == NULL) {
@ -114,13 +112,3 @@ double ComputeSMDTLSPHStrainRate::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}
/*
* deviator of a tensor
*/
Matrix3d ComputeSMDTLSPHStrainRate::Deviator(Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}

3
src/USER-SMD/compute_smd_tlsph_strain_rate.h
View File

@ -33,8 +33,6 @@ ComputeStyle(smd/tlsph_strain_rate,ComputeSMDTLSPHStrainRate)
#define LMP_COMPUTE_SMD_TLSPH_STRAIN_RATE_H
#include "compute.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
@ -45,7 +43,6 @@ class ComputeSMDTLSPHStrainRate : public Compute {
void init();
void compute_peratom();
double memory_usage();
Matrix3d Deviator(Matrix3d);
private:
int nmax;

21
src/USER-SMD/compute_smd_tlsph_stress.cpp
View File

@ -36,6 +36,17 @@
using namespace Eigen;
using namespace LAMMPS_NS;
/*
* deviator of a tensor
*/
static Matrix3d Deviator(Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHStress::ComputeSMDTLSPHStress(LAMMPS *lmp, int narg, char **arg) :
@ -119,13 +130,3 @@ double ComputeSMDTLSPHStress::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}
/*
* deviator of a tensor
*/
Matrix3d ComputeSMDTLSPHStress::Deviator(Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}

3
src/USER-SMD/compute_smd_tlsph_stress.h
View File

@ -33,8 +33,6 @@ ComputeStyle(smd/tlsph_stress, ComputeSMDTLSPHStress)
#define LMP_COMPUTE_SMD_TLSPH_STRESS_H
#include "compute.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
@ -45,7 +43,6 @@ class ComputeSMDTLSPHStress : public Compute {
void init();
void compute_peratom();
double memory_usage();
Matrix3d Deviator(Matrix3d);
private:
int nmax;

11
src/USER-SMD/compute_smd_ulsph_strain_rate.cpp
View File

@ -36,6 +36,7 @@
using namespace Eigen;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDULSPHStrainRate::ComputeSMDULSPHStrainRate(LAMMPS *lmp, int narg, char **arg) :
@ -120,13 +121,3 @@ double ComputeSMDULSPHStrainRate::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}
/*
* deviator of a tensor
*/
Matrix3d ComputeSMDULSPHStrainRate::Deviator(Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}

3
src/USER-SMD/compute_smd_ulsph_strain_rate.h
View File

@ -33,8 +33,6 @@ ComputeStyle(smd/ulsph_strain_rate,ComputeSMDULSPHStrainRate)
#define LMP_COMPUTE_SMD_ULSPH_STRAIN_RATE_H
#include "compute.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
@ -45,7 +43,6 @@ class ComputeSMDULSPHStrainRate : public Compute {
void init();
void compute_peratom();
double memory_usage();
Matrix3d Deviator(Matrix3d);
private:
int nmax;

19
src/USER-SMD/compute_smd_ulsph_stress.cpp
View File

@ -36,6 +36,16 @@
using namespace Eigen;
using namespace LAMMPS_NS;
/*
* deviator of a tensor
*/
static Matrix3d Deviator(Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}
/* ---------------------------------------------------------------------- */
ComputeSMDULSPHStress::ComputeSMDULSPHStress(LAMMPS *lmp, int narg, char **arg) :
@ -122,12 +132,3 @@ double ComputeSMDULSPHStress::memory_usage() {
return bytes;
}
/*
* deviator of a tensor
*/
Matrix3d ComputeSMDULSPHStress::Deviator(Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}

3
src/USER-SMD/compute_smd_ulsph_stress.h
View File

@ -33,8 +33,6 @@ ComputeStyle(smd/ulsph_stress, ComputeSMDULSPHStress)
#define LMP_COMPUTE_SMD_ULSPH_STRESS_H
#include "compute.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
@ -45,7 +43,6 @@ class ComputeSMDULSPHStress : public Compute {
void init();
void compute_peratom();
double memory_usage();
Matrix3d Deviator(Matrix3d);
private:
int nmax;

8
src/USER-SMD/fix_smd_move_triangulated_surface.h
View File

@ -31,10 +31,8 @@ FixStyle(smd/move_tri_surf,FixSMDMoveTriSurf)
#ifndef LMP_FIX_SMD_INTEGRATE_TRIANGULAR_SURFACE_H
#define LMP_FIX_SMD_INTEGRATE_TRIANGULAR_SURFACE_H
#include "fix.h"
#include <Eigen/Eigen>
using namespace Eigen;
#include "fix.h"
namespace LAMMPS_NS {
@ -53,9 +51,9 @@ protected:
double dtv;
bool linearFlag, wiggleFlag, rotateFlag;
double vx, vy, vz;
Vector3d rotation_axis, origin;
Eigen::Vector3d rotation_axis, origin;
double rotation_period;
Matrix3d u_cross, uxu;
Eigen::Matrix3d u_cross, uxu;
double wiggle_travel, wiggle_max_travel, wiggle_direction;
};

6
src/USER-SMD/fix_smd_tlsph_reference_configuration.cpp
View File

@ -202,7 +202,6 @@ void FixSMD_TLSPH_ReferenceConfiguration::pre_exchange() {
void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) {
int i, j, ii, jj, n, inum, jnum;
int *ilist, *jlist, *numneigh, **firstneigh;
int itype, jtype;
double r, h, wf, wfd;
Vector3d dx;
@ -219,7 +218,6 @@ void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) {
double **x0 = atom->x;
double *radius = atom->radius;
int *type = atom->type;
int *mask = atom->mask;
tagint *tag = atom->tag;
NeighList *list = pair->list;
@ -234,14 +232,12 @@ void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) {
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
if (INSERT_PREDEFINED_CRACKS) {
if (!crack_exclude(i, j))
@ -284,7 +280,6 @@ void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) {
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
@ -296,7 +291,6 @@ void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) {
dx(1) = x0[i][1] - x0[j][1];
dx(2) = x0[i][2] - x0[j][2];
r = dx.norm();
jtype = type[j];
h = radius[i] + radius[j];
if (INSERT_PREDEFINED_CRACKS) {

12
src/USER-SMD/fix_smd_wall_surface.cpp
View File

@ -31,8 +31,8 @@
#include "error.h"
#include <Eigen/Eigen>
#include <stdio.h>
#include <iostream>
#include "atom_vec.h"
#include <string.h>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -60,7 +60,7 @@ FixSMDWallSurface::FixSMDWallSurface(LAMMPS *lmp, int narg, char **arg) :
if (narg != 6)
error->all(FLERR, "Illegal number of arguments for fix smd/wall_surface");
filename.assign(arg[3]);
filename = strdup(arg[3]);
wall_particle_type = force->inumeric(FLERR, arg[4]);
wall_molecule_id = force->inumeric(FLERR, arg[5]);
if (wall_molecule_id < 65535) {
@ -69,7 +69,7 @@ FixSMDWallSurface::FixSMDWallSurface(LAMMPS *lmp, int narg, char **arg) :
if (comm->me == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
printf("fix smd/wall_surface reads trianglulated surface from file: %s\n", filename.c_str());
printf("fix smd/wall_surface reads trianglulated surface from file: %s\n", filename);
printf("fix smd/wall_surface has particle type %d \n", wall_particle_type);
printf("fix smd/wall_surface has molecule id %d \n", wall_molecule_id);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n");
@ -79,6 +79,8 @@ FixSMDWallSurface::FixSMDWallSurface(LAMMPS *lmp, int narg, char **arg) :
/* ---------------------------------------------------------------------- */
FixSMDWallSurface::~FixSMDWallSurface() {
free(filename);
filename = NULL;
// unregister this fix so atom class doesn't invoke it any more
//atom->delete_callback(id, 0);
@ -237,10 +239,10 @@ void FixSMDWallSurface::read_triangles(int pass) {
vert = new Vector3d[3];
Vector3d normal, center;
FILE *fp = fopen(filename.c_str(), "r");
FILE *fp = fopen(filename, "r");
if (fp == NULL) {
char str[128];
sprintf(str, "Cannot open file %s", filename.c_str());
sprintf(str, "Cannot open file %s", filename);
error->one(FLERR, str);
}

4
src/USER-SMD/fix_smd_wall_surface.h
View File

@ -21,8 +21,6 @@ FixStyle(smd/wall_surface,FixSMDWallSurface)
#define LMP_FIX_SMD_WALL_SURFACE_H
#include "fix.h"
#include <iostream>
using namespace std;
namespace LAMMPS_NS {
@ -42,7 +40,7 @@ public:
private:
int first; // flag for first time initialization
double sublo[3], subhi[3]; // epsilon-extended proc sub-box for adding atoms;
std::string filename;
char *filename;
int wall_particle_type;
int wall_molecule_id;
};

7
src/USER-SMD/pair_smd_tlsph.cpp
View File

@ -2102,7 +2102,6 @@ void PairTlsph::ComputeDamage(const int i, const Matrix3d strain, const Matrix3d
double *damage = atom->damage;
int *type = atom->type;
int itype = type[i];
bool damage_flag = false;
double jc_failure_strain;
//double damage_gap, damage_rate;
Matrix3d eye, stress_deviator;
@ -2116,16 +2115,15 @@ void PairTlsph::ComputeDamage(const int i, const Matrix3d strain, const Matrix3d
/*
* maximum stress failure criterion:
*/
damage_flag = IsotropicMaxStressDamage(stress, Lookup[FAILURE_MAX_PRINCIPAL_STRESS_THRESHOLD][itype]);
IsotropicMaxStressDamage(stress, Lookup[FAILURE_MAX_PRINCIPAL_STRESS_THRESHOLD][itype]);
} else if (failureModel[itype].failure_max_principal_strain) {
error->one(FLERR, "not yet implemented");
/*
* maximum strain failure criterion:
*/
damage_flag = IsotropicMaxStrainDamage(strain, Lookup[FAILURE_MAX_PRINCIPAL_STRAIN_THRESHOLD][itype]);
IsotropicMaxStrainDamage(strain, Lookup[FAILURE_MAX_PRINCIPAL_STRAIN_THRESHOLD][itype]);
} else if (failureModel[itype].failure_max_plastic_strain) {
if (eff_plastic_strain[i] >= Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype]) {
damage_flag = true;
damage[i] = 1.0;
//double damage_gap = 0.5 * Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype];
//damage[i] = (eff_plastic_strain[i] - Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype]) / damage_gap;
@ -2142,7 +2140,6 @@ void PairTlsph::ComputeDamage(const int i, const Matrix3d strain, const Matrix3d
//printf("JC failure strain is: %f\n", jc_failure_strain);
if (eff_plastic_strain[i] >= jc_failure_strain) {
damage_flag = true;
double damage_rate = Lookup[SIGNAL_VELOCITY][itype] / (100.0 * radius[i]);
damage[i] += damage_rate * update->dt;
//damage[i] = 1.0;

25
src/USER-SMD/pair_smd_tlsph.h
View File

@ -33,12 +33,7 @@ PairStyle(smd/tlsph,PairTlsph)
#include "pair.h"
#include <Eigen/Eigen>
#include <Eigen/Dense>
#include <Eigen/SVD>
#include <map>
using namespace std;
using namespace Eigen;
namespace LAMMPS_NS {
class PairTlsph: public Pair {
@ -67,13 +62,13 @@ public:
void PreCompute();
void ComputeForces(int eflag, int vflag);
void effective_longitudinal_modulus(const int itype, const double dt, const double d_iso, const double p_rate,
const Matrix3d d_dev, const Matrix3d sigma_dev_rate, const double damage, double &K_eff, double &mu_eff, double &M_eff);
const Eigen::Matrix3d d_dev, const Eigen::Matrix3d sigma_dev_rate, const double damage, double &K_eff, double &mu_eff, double &M_eff);
void ComputePressure(const int i, const double rho, const double mass_specific_energy, const double vol_specific_energy,
const double pInitial, const double d_iso, double &pFinal, double &p_rate);
void ComputeStressDeviator(const int i, const Matrix3d sigmaInitial_dev, const Matrix3d d_dev, Matrix3d &sigmaFinal_dev,
Matrix3d &sigma_dev_rate, double &plastic_strain_increment);
void ComputeDamage(const int i, const Matrix3d strain, const Matrix3d sigmaFinal, Matrix3d &sigma_damaged);
void ComputeStressDeviator(const int i, const Eigen::Matrix3d sigmaInitial_dev, const Eigen::Matrix3d d_dev, Eigen::Matrix3d &sigmaFinal_dev,
Eigen::Matrix3d &sigma_dev_rate, double &plastic_strain_increment);
void ComputeDamage(const int i, const Eigen::Matrix3d strain, const Eigen::Matrix3d sigmaFinal, Eigen::Matrix3d &sigma_damaged);
protected:
@ -89,12 +84,12 @@ protected:
/*
* per atom arrays
*/
Matrix3d *K, *PK1, *Fdot, *Fincr;
Matrix3d *R; // rotation matrix
Matrix3d *FincrInv;
Matrix3d *D, *W; // strain rate and spin tensor
Vector3d *smoothVelDifference;
Matrix3d *CauchyStress;
Eigen::Matrix3d *K, *PK1, *Fdot, *Fincr;
Eigen::Matrix3d *R; // rotation matrix
Eigen::Matrix3d *FincrInv;
Eigen::Matrix3d *D, *W; // strain rate and spin tensor
Eigen::Vector3d *smoothVelDifference;
Eigen::Matrix3d *CauchyStress;
double *detF, *particle_dt;
double *hourglass_error;
int *numNeighsRefConfig;

5
src/USER-SMD/pair_smd_triangulated_surface.h
View File

@ -34,7 +34,6 @@ PairStyle(smd/tri_surface,PairTriSurf)
#include "pair.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
@ -50,8 +49,8 @@ class PairTriSurf : public Pair {
void init_style();
void init_list(int, class NeighList *);
virtual double memory_usage();
void PointTriangleDistance(const Vector3d P, const Vector3d TRI1, const Vector3d TRI2, const Vector3d TRI3,
Vector3d &CP, double &dist);
void PointTriangleDistance(const Eigen::Vector3d P, const Eigen::Vector3d TRI1, const Eigen::Vector3d TRI2, const Eigen::Vector3d TRI3,
Eigen::Vector3d &CP, double &dist);
double clamp(const double a, const double min, const double max);
void *extract(const char *, int &);

5
src/USER-SMD/pair_smd_ulsph.cpp
View File

@ -222,7 +222,7 @@ void PairULSPH::PreCompute() {
int *ilist, *jlist, *numneigh;
int **firstneigh;
int nlocal = atom->nlocal;
int inum, jnum, ii, jj, i, itype, jtype, j, idim;
int inum, jnum, ii, jj, i, itype, j, idim;
double wfd, h, irad, r, rSq, wf, ivol, jvol;
Vector3d dx, dv, g, du;
Matrix3d Ktmp, Ltmp, Ftmp, K3di, D;
@ -281,7 +281,6 @@ void PairULSPH::PreCompute() {
if (rSq < h * h) {
r = sqrt(rSq);
jtype = type[j];
jvol = vfrac[j];
// distance vectors in current and reference configuration, velocity difference
@ -649,7 +648,6 @@ void PairULSPH::AssembleStressTensor() {
double **tlsph_stress = atom->smd_stress;
double *e = atom->e;
int *type = atom->type;
double pFinal;
int i, itype;
int nlocal = atom->nlocal;
Matrix3d D, Ddev, W, V, sigma_diag;
@ -686,7 +684,6 @@ void PairULSPH::AssembleStressTensor() {
error->one(FLERR, "unknown EOS.");
break;
case NONE:
pFinal = 0.0;
c0[i] = 1.0;
break;
case EOS_TAIT:

16
src/USER-SMD/pair_smd_ulsph.h
View File

@ -33,11 +33,7 @@ PairStyle(smd/ulsph,PairULSPH)
#include "pair.h"
#include <Eigen/Eigen>
#include <Eigen/Dense>
#include <Eigen/SVD>
using namespace Eigen;
using namespace std;
namespace LAMMPS_NS {
class PairULSPH: public Pair {
@ -57,10 +53,10 @@ public:
void *extract(const char *, int &);
void PreCompute();
void PreCompute_DensitySummation();
double effective_shear_modulus(const Matrix3d d_dev, const Matrix3d deltaStressDev, const double dt, const int itype);
double effective_shear_modulus(const Eigen::Matrix3d d_dev, const Eigen::Matrix3d deltaStressDev, const double dt, const int itype);
Vector3d ComputeHourglassForce(const int i, const int itype, const int j, const int jtype, const Vector3d dv,
const Vector3d xij, const Vector3d g, const double c_ij, const double mu_ij, const double rho_ij);
Eigen::Vector3d ComputeHourglassForce(const int i, const int itype, const int j, const int jtype, const Eigen::Vector3d dv,
const Eigen::Vector3d xij, const Eigen::Vector3d g, const double c_ij, const double mu_ij, const double rho_ij);
protected:
@ -78,10 +74,10 @@ protected:
int nmax; // max number of atoms on this proc
int *numNeighs;
Matrix3d *K;
Eigen::Matrix3d *K;
double *shepardWeight, *c0, *rho;
Vector3d *smoothVel;
Matrix3d *stressTensor, *L, *F;
Eigen::Vector3d *smoothVel;
Eigen::Matrix3d *stressTensor, *L, *F;
double dtCFL;