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

This commit is contained in:
sjplimp 2014-05-02 21:38:42 +00:00
parent 1fb3407215
commit ddcab1268e
29 changed files with 409 additions and 5010 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

427
src/USER-REAXC/fix_qeq_reax.cpp
View File

@ -14,6 +14,8 @@
/* ----------------------------------------------------------------------
Contributing author: Hasan Metin Aktulga, Purdue University
(now at Lawrence Berkeley National Laboratory, hmaktulga@lbl.gov)
Hybrid and sub-group capabilities: Ray Shan (Sandia)
------------------------------------------------------------------------- */
#include "math.h"
@ -30,6 +32,7 @@
#include "neigh_request.h"
#include "update.h"
#include "force.h"
#include "group.h"
#include "pair.h"
#include "respa.h"
#include "memory.h"
@ -68,8 +71,6 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) :
if (narg != 8) error->all(FLERR,"Illegal fix qeq/reax command");
nevery = force->inumeric(FLERR,arg[3]);
if (nevery <= 0) error->all(FLERR,"Illegal fix qeq/reax command");
swa = force->numeric(FLERR,arg[4]);
swb = force->numeric(FLERR,arg[5]);
tolerance = force->numeric(FLERR,arg[6]);
@ -103,15 +104,6 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) :
H.jlist = NULL;
H.val = NULL;
// GMRES
//g = NULL;
//y = NULL;
//hstr = NULL;
//v = NULL;
//h = NULL;
//hc = NULL;
//hs = NULL;
// perform initial allocation of atom-based arrays
// register with Atom class
@ -124,6 +116,7 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) :
reaxc = NULL;
reaxc = (PairReaxC *) force->pair_match("reax/c",1);
}
/* ---------------------------------------------------------------------- */
@ -262,7 +255,8 @@ void FixQEqReax::reallocate_storage()
void FixQEqReax::allocate_matrix()
{
int i,ii;
int i,ii,inum,m;
int *ilist, *numneigh;
int mincap;
double safezone;
@ -280,10 +274,20 @@ void FixQEqReax::allocate_matrix()
// determine the total space for the H matrix
int m = 0;
for( ii = 0; ii < list->inum; ii++ ) {
i = list->ilist[ii];
m += list->numneigh[i];
if (reaxc) {
inum = reaxc->list->inum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
} else {
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
}
m = 0;
for( ii = 0; ii < inum; ii++ ) {
i = ilist[ii];
m += numneigh[i];
}
m_cap = MAX( (int)(m * safezone), mincap * MIN_NBRS );
@ -322,15 +326,24 @@ void FixQEqReax::init()
if (!force->pair_match("reax/c",1))
error->all(FLERR,"Must use pair_style reax/c with fix qeq/reax");
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/reax group has no atoms");
/*
if (reaxc)
if (ngroup != reaxc->ngroup)
error->all(FLERR,"Fix qeq/reax group and pair reax/c have "
"different numbers of atoms");
*/
// need a half neighbor list w/ Newton off and ghost neighbors
// make it occasional if QeQ not performed every timestep
// built whenever re-neighboring occurs
int irequest = neighbor->request(this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->newton = 2;
neighbor->requests[irequest]->ghost = 1;
if (nevery > 1) neighbor->requests[irequest]->occasional = 1;
init_shielding();
init_taper();
@ -427,8 +440,14 @@ void FixQEqReax::min_setup_pre_force(int vflag)
void FixQEqReax::init_storage()
{
N = atom->nlocal + atom->nghost;
for( int i = 0; i < N; i++ ) {
int NN;
if (reaxc)
NN = reaxc->list->inum + reaxc->list->gnum;
else
NN = list->inum + list->gnum;
for( int i = 0; i < NN; i++ ) {
Hdia_inv[i] = 1. / eta[atom->type[i]];
b_s[i] = -chi[atom->type[i]];
b_t[i] = -1.0;
@ -449,15 +468,17 @@ void FixQEqReax::pre_force(int vflag)
n = atom->nlocal;
N = atom->nlocal + atom->nghost;
// grow arrays if necessary
// need to be atom->nmax in length
if( atom->nmax > nmax ) reallocate_storage();
if( n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE )
reallocate_matrix();
init_matvec();
matvecs = CG(b_s, s); // CG on s - parallel
matvecs += CG(b_t, t); // CG on t - parallel
matvecs = CG(b_s, s); // CG on s - parallel
matvecs += CG(b_t, t); // CG on t - parallel
calculate_Q();
if( comm->me == 0 ) {
@ -487,23 +508,38 @@ void FixQEqReax::init_matvec()
/* fill-in H matrix */
compute_H();
for( int i = 0; i < n; ++i ) {
/* init pre-conditioner for H and init solution vectors */
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -chi[ atom->type[i] ];
b_t[i] = -1.0;
int nn, ii, i;
int *ilist;
/* linear extrapolation for s & t from previous solutions */
//s[i] = 2 * s_hist[i][0] - s_hist[i][1];
//t[i] = 2 * t_hist[i][0] - t_hist[i][1];
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
/* quadratic extrapolation for s & t from previous solutions */
//s[i] = s_hist[i][2] + 3 * ( s_hist[i][0] - s_hist[i][1] );
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] );
for( ii = 0; ii < nn; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
/* init pre-conditioner for H and init solution vectors */
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -chi[ atom->type[i] ];
b_t[i] = -1.0;
/* cubic extrapolation for s & t from previous solutions */
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
//t[i] = 4*(t_hist[i][0]+t_hist[i][2])-(6*t_hist[i][1]+t_hist[i][3]);
/* linear extrapolation for s & t from previous solutions */
//s[i] = 2 * s_hist[i][0] - s_hist[i][1];
//t[i] = 2 * t_hist[i][0] - t_hist[i][1];
/* quadratic extrapolation for s & t from previous solutions */
//s[i] = s_hist[i][2] + 3 * ( s_hist[i][0] - s_hist[i][1] );
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] );
/* cubic extrapolation for s & t from previous solutions */
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
//t[i] = 4*(t_hist[i][0]+t_hist[i][2])-(6*t_hist[i][1]+t_hist[i][3]);
}
}
pack_flag = 2;
@ -517,60 +553,67 @@ void FixQEqReax::init_matvec()
void FixQEqReax::compute_H()
{
int inum, jnum, *ilist, *jlist, *numneigh, **firstneigh;
int i, j, ii, jj, flag;
int i, j, ii, jj, temp, newnbr, flag;
double **x, SMALL = 0.0001;
double dx, dy, dz, r_sqr;
int *type = atom->type;
tagint *tag = atom->tag;
x = atom->x;
int *mask = atom->mask;
if (nevery > 1) neighbor->build_one(list->index);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
if (reaxc) {
inum = reaxc->list->inum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
// fill in the H matrix
m_fill = 0;
r_sqr = 0;
for( ii = 0; ii < inum; ii++ ) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
H.firstnbr[i] = m_fill;
if (mask[i] & groupbit) {
jlist = firstneigh[i];
jnum = numneigh[i];
H.firstnbr[i] = m_fill;
for( jj = 0; jj < jnum; jj++ ) {
j = jlist[jj];
for( jj = 0; jj < jnum; jj++ ) {
j = jlist[jj];
dx = x[j][0] - x[i][0];
dy = x[j][1] - x[i][1];
dz = x[j][2] - x[i][2];
r_sqr = SQR(dx) + SQR(dy) + SQR(dz);
dx = x[j][0] - x[i][0];
dy = x[j][1] - x[i][1];
dz = x[j][2] - x[i][2];
r_sqr = SQR(dx) + SQR(dy) + SQR(dz);
flag = 0;
if (r_sqr <= SQR(swb)) {
if (j < n) flag = 1;
else if (tag[i] < tag[j]) flag = 1;
else if (tag[i] == tag[j]) {
if (dz > SMALL) flag = 1;
else if (fabs(dz) < SMALL) {
if (dy > SMALL) flag = 1;
else if (fabs(dy) < SMALL && dx > SMALL)
flag = 1;
}
flag = 0;
if (r_sqr <= SQR(swb)) {
if (j < n) flag = 1;
else if (tag[i] < tag[j]) flag = 1;
else if (tag[i] == tag[j]) {
if (dz > SMALL) flag = 1;
else if (fabs(dz) < SMALL) {
if (dy > SMALL) flag = 1;
else if (fabs(dy) < SMALL && dx > SMALL)
flag = 1;
}
}
}
if( flag ) {
H.jlist[m_fill] = j;
H.val[m_fill] = calculate_H( sqrt(r_sqr), shld[type[i]][type[j]] );
m_fill++;
}
}
if( flag ) {
H.jlist[m_fill] = j;
H.val[m_fill] = calculate_H( sqrt(r_sqr), shld[type[i]][type[j]] );
m_fill++;
}
H.numnbrs[i] = m_fill - H.firstnbr[i];
}
H.numnbrs[i] = m_fill - H.firstnbr[i];
}
if (m_fill >= H.m) {
@ -608,7 +651,17 @@ int FixQEqReax::CG( double *b, double *x )
{
int i, j, imax;
double tmp, alpha, beta, b_norm;
double sig_old, sig_new;
double sig_old, sig_new, sig0;
int nn, jj;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
imax = 200;
@ -616,35 +669,44 @@ int FixQEqReax::CG( double *b, double *x )
sparse_matvec( &H, x, q );
comm->reverse_comm_fix( this ); //Coll_Vector( q );
vector_sum( r , 1., b, -1., q, n );
for( j = 0; j < n; ++j )
d[j] = r[j] * Hdia_inv[j]; //pre-condition
vector_sum( r , 1., b, -1., q, nn );
b_norm = parallel_norm( b, n );
sig_new = parallel_dot( r, d, n );
for( jj = 0; jj < nn; ++jj ) {
j = ilist[jj];
if (atom->mask[j] & groupbit)
d[j] = r[j] * Hdia_inv[j]; //pre-condition
}
int ttype = 1;
b_norm = parallel_norm( b, nn );
sig_new = parallel_dot( r, d, nn);
sig0 = sig_new;
for( i = 1; i < imax && sqrt(sig_new) / b_norm > tolerance; ++i ) {
comm->forward_comm_fix(this); //Dist_vector( d );
sparse_matvec( &H, d, q );
comm->reverse_comm_fix(this); //Coll_vector( q );
tmp = parallel_dot( d, q, n );
ttype = 2;
tmp = parallel_dot( d, q, nn);
alpha = sig_new / tmp;
// comm->me, i, parallel_norm( d, n ), parallel_norm( q, n ), tmp );
vector_add( x, alpha, d, n );
vector_add( r, -alpha, q, n );
vector_add( x, alpha, d, nn );
vector_add( r, -alpha, q, nn );
// pre-conditioning
for( j = 0; j < n; ++j )
p[j] = r[j] * Hdia_inv[j];
for( jj = 0; jj < nn; ++jj ) {
j = ilist[jj];
if (atom->mask[j] & groupbit)
p[j] = r[j] * Hdia_inv[j];
}
sig_old = sig_new;
sig_new = parallel_dot( r, p, n );
sig_new = parallel_dot( r, p, nn);
beta = sig_new / sig_old;
vector_sum( d, 1., p, beta, d, n );
vector_sum( d, 1., p, beta, d, nn );
}
if (i >= imax && comm->me == 0) {
@ -663,19 +725,42 @@ int FixQEqReax::CG( double *b, double *x )
void FixQEqReax::sparse_matvec( sparse_matrix *A, double *x, double *b )
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
for( i = 0; i < n; ++i )
b[i] = eta[ atom->type[i] ] * x[i];
for( i = n; i < N; ++i )
b[i] = 0;
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
}
for( i = 0; i < n; ++i ) {
for( itr_j=A->firstnbr[i]; itr_j<A->firstnbr[i]+A->numnbrs[i]; itr_j++) {
j = A->jlist[itr_j];
b[i] += A->val[itr_j] * x[j];
b[j] += A->val[itr_j] * x[i];
for( ii = 0; ii < nn; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
b[i] = eta[ atom->type[i] ] * x[i];
}
for( ii = nn; ii < NN; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
b[i] = 0;
}
for( ii = 0; ii < nn; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
for( itr_j=A->firstnbr[i]; itr_j<A->firstnbr[i]+A->numnbrs[i]; itr_j++) {
j = A->jlist[itr_j];
b[i] += A->val[itr_j] * x[j];
b[j] += A->val[itr_j] * x[i];
}
}
}
}
/* ---------------------------------------------------------------------- */
@ -686,20 +771,34 @@ void FixQEqReax::calculate_Q()
double u, s_sum, t_sum;
double *q = atom->q;
s_sum = parallel_vector_acc( s, n );
t_sum = parallel_vector_acc( t, n);
int nn, ii;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
s_sum = parallel_vector_acc( s, nn );
t_sum = parallel_vector_acc( t, nn);
u = s_sum / t_sum;
for( i = 0; i < n; ++i ) {
q[i] = s[i] - u * t[i];
for( ii = 0; ii < nn; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
q[i] = s[i] - u * t[i];
/* backup s & t */
for( k = 4; k > 0; --k ) {
s_hist[i][k] = s_hist[i][k-1];
t_hist[i][k] = t_hist[i][k-1];
/* backup s & t */
for( k = 4; k > 0; --k ) {
s_hist[i][k] = s_hist[i][k-1];
t_hist[i][k] = t_hist[i][k-1];
}
s_hist[i][0] = s[i];
t_hist[i][0] = t[i];
}
s_hist[i][0] = s[i];
t_hist[i][0] = t[i];
}
pack_flag = 4;
@ -825,9 +924,21 @@ double FixQEqReax::parallel_norm( double *v, int n )
int i;
double my_sum, norm_sqr;
my_sum = 0;
for( i = 0; i < n; ++i )
my_sum += SQR( v[i] );
int ii;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
my_sum = 0.0;
norm_sqr = 0.0;
for( ii = 0; ii < n; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
my_sum += SQR( v[i] );
}
MPI_Allreduce( &my_sum, &norm_sqr, 1, MPI_DOUBLE, MPI_SUM, world );
@ -836,16 +947,27 @@ double FixQEqReax::parallel_norm( double *v, int n )
/* ---------------------------------------------------------------------- */
double FixQEqReax::parallel_dot( double *v1, double *v2, int n )
double FixQEqReax::parallel_dot( double *v1, double *v2, int n)
{
int i;
double my_dot, res;
my_dot = 0;
res = 0;
for( i = 0; i < n; ++i )
my_dot += v1[i] * v2[i];
int ii;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
my_dot = 0.0;
res = 0.0;
for( ii = 0; ii < n; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
my_dot += v1[i] * v2[i];
}
MPI_Allreduce( &my_dot, &res, 1, MPI_DOUBLE, MPI_SUM, world );
return res;
@ -858,9 +980,21 @@ double FixQEqReax::parallel_vector_acc( double *v, int n )
int i;
double my_acc, res;
my_acc = 0;
for( i = 0; i < n; ++i )
my_acc += v[i];
int ii;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
my_acc = 0.0;
res = 0.0;
for( ii = 0; ii < n; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit)
my_acc += v[i];
}
MPI_Allreduce( &my_acc, &res, 1, MPI_DOUBLE, MPI_SUM, world );
@ -869,49 +1003,40 @@ double FixQEqReax::parallel_vector_acc( double *v, int n )
/* ---------------------------------------------------------------------- */
double FixQEqReax::norm( double* v1, int k )
{
double ret = 0;
for( --k; k>=0; --k )
ret += ( v1[k] * v1[k] );
return sqrt( ret );
}
/* ---------------------------------------------------------------------- */
void FixQEqReax::vector_sum( double* dest, double c, double* v,
double d, double* y, int k )
{
for( --k; k>=0; --k )
dest[k] = c * v[k] + d * y[k];
}
int kk;
int *ilist;
/* ---------------------------------------------------------------------- */
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
void FixQEqReax::vector_scale( double* dest, double c, double* v, int k )
{
for( --k; k>=0; --k )
dest[k] = c * v[k];
}
/* ---------------------------------------------------------------------- */
double FixQEqReax::dot( double* v1, double* v2, int k )
{
double ret = 0;
for( --k; k>=0; --k )
ret += v1[k] * v2[k];
return ret;
for( --k; k>=0; --k ) {
kk = ilist[k];
if (atom->mask[kk] & groupbit)
dest[kk] = c * v[kk] + d * y[kk];
}
}
/* ---------------------------------------------------------------------- */
void FixQEqReax::vector_add( double* dest, double c, double* v, int k )
{
for( --k; k>=0; --k )
dest[k] += c * v[k];
int kk;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
for( --k; k>=0; --k ) {
kk = ilist[k];
if (atom->mask[kk] & groupbit)
dest[kk] += c * v[kk];
}
}

5
src/USER-REAXC/fix_qeq_reax.h
View File

@ -70,6 +70,8 @@ class FixQEqReax : public Fix {
double *chi,*eta,*gamma; // qeq parameters
double **shld;
bigint ngroup;
// fictitious charges
double *s, *t;
@ -132,10 +134,7 @@ class FixQEqReax : public Fix {
double parallel_dot( double*, double*, int );
double parallel_vector_acc( double*, int );
double norm(double*,int);
void vector_sum(double*,double,double*,double,double*,int);
void vector_scale(double*,double,double*,int);
double dot(double*,double*,int);
void vector_add(double*, double, double*,int);
};

1
src/USER-REAXC/fix_reaxc_species.cpp
View File

@ -53,6 +53,7 @@ FixReaxCSpecies::FixReaxCSpecies(LAMMPS *lmp, int narg, char **arg) :
vector_flag = 1;
size_vector = 2;
extvector = 0;
peratom_flag = 1;
size_peratom_cols = 0;

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

@ -17,6 +17,7 @@
Per-atom energy/virial added by Ray Shan (Sandia)
Fix reax/c/bonds and fix reax/c/species for pair_style reax/c added by
Ray Shan (Sandia)
Hybrid and hybrid/overlay compatibility added by Ray Shan (Sandia)
------------------------------------------------------------------------- */
#include "pair_reax_c.h"
@ -283,11 +284,12 @@ void PairReaxC::coeff( int nargs, char **args )
// read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL
int itmp;
int itmp = 0;
int nreax_types = system->reax_param.num_atom_types;
for (int i = 3; i < nargs; i++) {
if (strcmp(args[i],"NULL") == 0) {
map[i-2] = -1;
itmp ++;
continue;
}
}
@ -295,7 +297,6 @@ void PairReaxC::coeff( int nargs, char **args )
int n = atom->ntypes;
// pair_coeff element map
itmp = 0;
for (int i = 3; i < nargs; i++)
for (int j = 0; j < nreax_types; j++)
if (strcasecmp(args[i],system->reax_param.sbp[j].name) == 0) {
@ -307,14 +308,22 @@ void PairReaxC::coeff( int nargs, char **args )
if (itmp != n)
error->all(FLERR,"Non-existent ReaxFF type");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++) {
setflag[i][j] = 1;
count++;
}
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ---------------------------------------------------------------------- */
@ -367,6 +376,7 @@ void PairReaxC::init_style( )
delete [] fixarg;
fix_reax = (FixReaxC *) modify->fix[modify->nfix-1];
}
}
/* ---------------------------------------------------------------------- */
@ -429,12 +439,23 @@ void PairReaxC::setup( )
ReAllocate( system, control, data, workspace, &lists, mpi_data );
}
ngroup = 0;
int ngroup_sum = 0;
for (int i = 0; i < list->inum; i++) {
ngroup ++;
}
MPI_Allreduce( &ngroup, &ngroup_sum, 1, MPI_INT, MPI_SUM, world );
ngroup = ngroup_sum;
}
/* ---------------------------------------------------------------------- */
double PairReaxC::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
cutghost[i][j] = cutghost[j][i] = cutmax;
return cutmax;
}
@ -485,7 +506,7 @@ void PairReaxC::compute(int eflag, int vflag)
// forces
Compute_Forces(system,control,data,workspace,&lists,out_control,mpi_data);
read_reax_forces();
read_reax_forces(vflag);
for(int k = 0; k < system->N; ++k) {
num_bonds[k] = system->my_atoms[k].num_bonds;
@ -609,12 +630,12 @@ void PairReaxC::set_far_nbr( far_neighbor_data *fdest,
int PairReaxC::estimate_reax_lists()
{
int itr_i, itr_j, i, j;
int itr_i, itr_j, itr_g, i, j, g;
int nlocal, nghost, num_nbrs, num_marked;
int *ilist, *jlist, *numneigh, **firstneigh, *marked;
double d_sqr;
rvec dvec;
double *dist, **x;
double d_sqr, g_d_sqr;
rvec dvec, g_dvec;
double **x;
reax_list *far_nbrs;
far_neighbor_data *far_list;
@ -634,11 +655,12 @@ int PairReaxC::estimate_reax_lists()
num_nbrs = 0;
num_marked = 0;
marked = (int*) calloc( system->N, sizeof(int) );
dist = (double*) calloc( system->N, sizeof(double) );
int numall = list->inum + list->gnum;
int inum = list->inum;
int gnum = list->gnum;
int numall = inum + gnum;
for( itr_i = 0; itr_i < numall; ++itr_i ){
for( itr_i = 0; itr_i < inum+gnum; ++itr_i ){
i = ilist[itr_i];
marked[i] = 1;
++num_marked;
@ -655,7 +677,6 @@ int PairReaxC::estimate_reax_lists()
}
free( marked );
free( dist );
return static_cast<int> (MAX( num_nbrs*safezone, mincap*MIN_NBRS ));
}
@ -664,12 +685,12 @@ int PairReaxC::estimate_reax_lists()
int PairReaxC::write_reax_lists()
{
int itr_i, itr_j, i, j;
int itr_i, itr_j, itr_g, i, j, g, flag;
int nlocal, nghost, num_nbrs;
int *ilist, *jlist, *numneigh, **firstneigh, *marked;
double d_sqr;
rvec dvec;
double *dist, **x;
int *ilist, *jlist, *numneigh, **firstneigh;
double d_sqr, g_d, g_d_sqr;
rvec dvec, g_dvec;
double *dist, **x, SMALL = 0.0001;
reax_list *far_nbrs;
far_neighbor_data *far_list;
@ -684,14 +705,14 @@ int PairReaxC::write_reax_lists()
far_list = far_nbrs->select.far_nbr_list;
num_nbrs = 0;
marked = (int*) calloc( system->N, sizeof(int) );
dist = (double*) calloc( system->N, sizeof(double) );
int numall = list->inum + list->gnum;
int inum = list->inum;
int gnum = list->gnum;
int numall = inum + gnum;
for( itr_i = 0; itr_i < numall; ++itr_i ){
for( itr_i = 0; itr_i < inum+gnum; ++itr_i ){
i = ilist[itr_i];
marked[i] = 1;
jlist = firstneigh[i];
Set_Start_Index( i, num_nbrs, far_nbrs );
@ -709,7 +730,6 @@ int PairReaxC::write_reax_lists()
Set_End_Index( i, num_nbrs, far_nbrs );
}
free( marked );
free( dist );
return num_nbrs;
@ -717,17 +737,18 @@ int PairReaxC::write_reax_lists()
/* ---------------------------------------------------------------------- */
void PairReaxC::read_reax_forces()
void PairReaxC::read_reax_forces(int vflag)
{
for( int i = 0; i < system->N; ++i ) {
system->my_atoms[i].f[0] = workspace->f[i][0];
system->my_atoms[i].f[1] = workspace->f[i][1];
system->my_atoms[i].f[2] = workspace->f[i][2];
atom->f[i][0] = -workspace->f[i][0];
atom->f[i][1] = -workspace->f[i][1];
atom->f[i][2] = -workspace->f[i][2];
atom->f[i][0] += -workspace->f[i][0];
atom->f[i][1] += -workspace->f[i][1];
atom->f[i][2] += -workspace->f[i][2];
}
}
/* ---------------------------------------------------------------------- */
@ -771,6 +792,8 @@ double PairReaxC::memory_usage()
bytes += 19.0 * system->total_cap * sizeof(real);
bytes += 3.0 * system->total_cap * sizeof(int);
double mem1 = bytes;
// From reaxc_lists
bytes += 2.0 * lists->n * sizeof(int);
bytes += lists->num_intrs * sizeof(three_body_interaction_data);
@ -790,8 +813,8 @@ double PairReaxC::memory_usage()
void PairReaxC::FindBond()
{
int i, j, pj, nj;
double bo_tmp, bo_cut, r_tmp;
int i, ii, j, pj, nj, jtmp, jj;
double bo_tmp, bo_cut, rij, rsq, r_tmp;
bond_data *bo_ij;
bo_cut = 0.10;
@ -815,5 +838,3 @@ void PairReaxC::FindBond()
}
}
}
/* ---------------------------------------------------------------------- */

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

@ -57,8 +57,12 @@ class PairReaxC : public Pair {
reax_list *lists;
mpi_datatypes *mpi_data;
bigint ngroup;
private:
double cutmax;
int nelements; // # of unique elements
char **elements; // names of unique elements
int *map;
class FixReaxC *fix_reax;
@ -68,17 +72,20 @@ class PairReaxC : public Pair {
int firstwarn;
void allocate();
void setup();
void create_compute();
void create_fix();
void write_reax_atoms();
void get_distance(rvec, rvec, double *, rvec *);
void set_far_nbr(far_neighbor_data *, int, double, rvec);
int estimate_reax_lists();
int write_reax_lists();
void read_reax_forces();
void setup();
void read_reax_forces(int);
int nmax;
void FindBond();
double memory_usage();
};
}

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

@ -25,20 +25,11 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "allocate.h"
#include "list.h"
#include "reset_tools.h"
#include "tool_box.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_allocate.h"
#include "reaxc_list.h"
#include "reaxc_reset_tools.h"
#include "reaxc_tool_box.h"
#include "reaxc_vector.h"
#endif
/* allocate space for my_atoms
important: we cannot know the exact number of atoms that will fall into a
@ -55,28 +46,9 @@ int PreAllocate_Space( reax_system *system, control_params *control,
system->local_cap = MAX( (int)(system->n * safezone), mincap );
system->total_cap = MAX( (int)(system->N * safezone), mincap );
#if defined(DEBUG)
fprintf( stderr, "p%d: local_cap=%d total_cap=%d\n",
system->my_rank, system->local_cap, system->total_cap );
#endif
system->my_atoms = (reax_atom*)
scalloc( system->total_cap, sizeof(reax_atom), "my_atoms", comm );
/* space for keeping restriction info, if any */
// not yet implemented in the parallel version!!!
// if( control->restrict_bonds ) {
// workspace->restricted = (int*)
// scalloc( system->local_cap, sizeof(int), "restricted_atoms", comm );
// workspace->restricted_list = (int**)
// scalloc( system->local_cap, sizeof(int*), "restricted_list", comm );
// for( i = 0; i < system->local_cap; ++i )
// workspace->restricted_list[i] = (int*)
// scalloc( MAX_RESTRICT, sizeof(int), "restricted_list[i]", comm );
// }
return SUCCESS;
}
@ -228,57 +200,12 @@ void DeAllocate_Workspace( control_params *control, storage *workspace )
sfree( workspace->p2, "p2" );
/* integrator */
// sfree( workspace->f_old );
sfree( workspace->v_const, "v_const" );
/*workspace->realloc.num_far = -1;
workspace->realloc.Htop = -1;
workspace->realloc.hbonds = -1;
workspace->realloc.bonds = -1;
workspace->realloc.num_3body = -1;
workspace->realloc.gcell_atoms = -1;*/
/* storage for analysis */
// if( control->molecular_analysis || control->diffusion_coef ) {
// sfree( workspace->mark, "mark" );
// sfree( workspace->old_mark, "old_mark" );
// }
// if( control->diffusion_coef )
// sfree( workspace->x_old, "x_old" );
/* force related storage */
sfree( workspace->f, "f" );
sfree( workspace->CdDelta, "CdDelta" );
#ifdef TEST_FORCES
sfree(workspace->dDelta, "dDelta" );
sfree( workspace->f_ele, "f_ele" );
sfree( workspace->f_vdw, "f_vdw" );
sfree( workspace->f_bo, "f_bo" );
sfree( workspace->f_be, "f_be" );
sfree( workspace->f_lp, "f_lp" );
sfree( workspace->f_ov, "f_ov" );
sfree( workspace->f_un, "f_un" );
sfree( workspace->f_ang, "f_ang" );
sfree( workspace->f_coa, "f_coa" );
sfree( workspace->f_pen, "f_pen" );
sfree( workspace->f_hb, "f_hb" );
sfree( workspace->f_tor, "f_tor" );
sfree( workspace->f_con, "f_con" );
sfree( workspace->f_tot, "f_tot" );
sfree( workspace->rcounts, "rcounts" );
sfree( workspace->displs, "displs" );
sfree( workspace->id_all, "id_all" );
sfree( workspace->f_all, "f_all" );
#endif
/* hbond storage */
//sfree( workspace->Hindex, "Hindex" );
//sfree( workspace->num_bonds );
//sfree( workspace->num_hbonds );
//sfree( workspace->hash, "hash" );
//sfree( workspace->rev_hash, "rev_hash" );
}
@ -328,7 +255,6 @@ int Allocate_Workspace( reax_system *system, control_params *control,
scalloc( total_cap, sizeof(int), "bond_mark", comm );
workspace->done_after = (int*)
scalloc( total_cap, sizeof(int), "done_after", comm );
// fprintf( stderr, "p%d: bond order storage\n", system->my_rank );
/* QEq storage */
workspace->Hdia_inv = (real*)
@ -371,61 +297,11 @@ int Allocate_Workspace( reax_system *system, control_params *control,
/* integrator storage */
workspace->v_const = (rvec*) smalloc( local_rvec, "v_const", comm );
/* storage for analysis */
// not yet implemented in the parallel version!!!
// if( control->molecular_analysis || control->diffusion_coef ) {
// workspace->mark = (int*) scalloc( local_cap, sizeof(int), "mark", comm );
// workspace->old_mark = (int*)
// scalloc( local_cap, sizeof(int), "old_mark", comm );
// }
// else
// workspace->mark = workspace->old_mark = NULL;
// if( control->diffusion_coef )
// workspace->x_old = (rvec*)
// scalloc( local_cap, sizeof(rvec), "x_old", comm );
// else workspace->x_old = NULL;
// /* force related storage */
workspace->f = (rvec*) scalloc( total_cap, sizeof(rvec), "f", comm );
workspace->CdDelta = (real*)
scalloc( total_cap, sizeof(real), "CdDelta", comm );
#ifdef TEST_FORCES
workspace->dDelta=(rvec*) smalloc( total_rvec, "dDelta", comm );
workspace->f_ele =(rvec*) smalloc( total_rvec, "f_ele", comm );
workspace->f_vdw =(rvec*) smalloc( total_rvec, "f_vdw", comm );
workspace->f_bo =(rvec*) smalloc( total_rvec, "f_bo", comm );
workspace->f_be =(rvec*) smalloc( total_rvec, "f_be", comm );
workspace->f_lp =(rvec*) smalloc( total_rvec, "f_lp", comm );
workspace->f_ov =(rvec*) smalloc( total_rvec, "f_ov", comm );
workspace->f_un =(rvec*) smalloc( total_rvec, "f_un", comm );
workspace->f_ang =(rvec*) smalloc( total_rvec, "f_ang", comm );
workspace->f_coa =(rvec*) smalloc( total_rvec, "f_coa", comm );
workspace->f_pen =(rvec*) smalloc( total_rvec, "f_pen", comm );
workspace->f_hb =(rvec*) smalloc( total_rvec, "f_hb", comm );
workspace->f_tor =(rvec*) smalloc( total_rvec, "f_tor", comm );
workspace->f_con =(rvec*) smalloc( total_rvec, "f_con", comm );
workspace->f_tot =(rvec*) smalloc( total_rvec, "f_tot", comm );
if( system->my_rank == MASTER_NODE ) {
workspace->rcounts = (int*)
smalloc( system->wsize*sizeof(int), "rcount", comm );
workspace->displs = (int*)
smalloc( system->wsize*sizeof(int), "displs", comm );
workspace->id_all = (int*)
smalloc( system->bigN*sizeof(int), "id_all", comm );
workspace->f_all = (rvec*)
smalloc( system->bigN*sizeof(rvec), "f_all", comm );
}
else{
workspace->rcounts = NULL;
workspace->displs = NULL;
workspace->id_all = NULL;
workspace->f_all = NULL;
}
#endif
return SUCCESS;
}
@ -477,11 +353,6 @@ int Reallocate_Matrix( sparse_matrix **H, int n, int m, char *name,
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "reallocating %s matrix, n = %d, m = %d\n", name, n, m );
fprintf( stderr, "memory allocated: %s = %dMB\n",
name, (int)(m * sizeof(sparse_matrix_entry) / (1024*1024)) );
#endif
return SUCCESS;
}
@ -497,9 +368,6 @@ int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds,
total_hbonds = 0;
for( i = 0; i < system->n; ++i )
if( (id = system->my_atoms[i].Hindex) >= 0 ) {
// commented out - already updated in validate_lists in forces.c
// system->my_atoms[i].num_hbonds = MAX(Num_Entries(id,hbonds)*SAFER_ZONE,
// MIN_HBONDS);
total_hbonds += system->my_atoms[i].num_hbonds;
}
total_hbonds = (int)(MAX( total_hbonds*saferzone, mincap*MIN_HBONDS ));
@ -526,8 +394,6 @@ int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
*est_3body = 0;
for( i = 0; i < system->N; ++i ){
*est_3body += SQR(system->my_atoms[i].num_bonds);
// commented out - already updated in validate_lists in forces.c
// system->my_atoms[i].num_bonds = MAX( Num_Entries(i,bonds)*2, MIN_BONDS );
*total_bonds += system->my_atoms[i].num_bonds;
}
*total_bonds = (int)(MAX( *total_bonds * safezone, mincap*MIN_BONDS ));
@ -567,13 +433,6 @@ int Estimate_GCell_Population( reax_system* system, MPI_Comm comm )
else if( c[d] < g->native_str[d] )
c[d] = g->native_str[d];
}
#if defined(DEBUG)
fprintf( stderr, "p%d bin_my_atoms: l:%d - atom%d @ %.5f %.5f %.5f" \
"--> cell: %d %d %d\n",
system->my_rank, l, atoms[l].orig_id,
atoms[l].x[0], atoms[l].x[1], atoms[l].x[2],
c[0], c[1], c[2] );
#endif
gc = &( g->cells[c[0]][c[1]][c[2]] );
gc->top++;
}
@ -585,18 +444,10 @@ int Estimate_GCell_Population( reax_system* system, MPI_Comm comm )
gc = &(g->cells[i][j][k]);
if( max_atoms < gc->top )
max_atoms = gc->top;
#if defined(DEBUG)
fprintf( stderr, "p%d gc[%d,%d,%d]->top=%d\n",
system->my_rank, i, j, k, gc->top );
#endif
}
my_max = (int)(MAX(max_atoms*SAFE_ZONE, MIN_GCELL_POPL));
MPI_Allreduce( &my_max, &all_max, 1, MPI_INT, MPI_MAX, comm );
#if defined(DEBUG)
fprintf( stderr, "p%d max_atoms=%d, my_max=%d, all_max=%d\n",
system->my_rank, max_atoms, my_max, all_max );
#endif
return all_max;
}
@ -650,18 +501,6 @@ void Allocate_Grid( reax_system *system, MPI_Comm comm )
for( k = g->native_str[2]; k < g->native_end[2]; ++k )
g->cells[i][j][k].atoms = (int*) scalloc( g->max_atoms, sizeof(int),
"g:atoms", comm );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d-allocated %dx%dx%d grid: nbrs=%d atoms=%d space=%dMB\n",
system->my_rank, g->ncells[0], g->ncells[1], g->ncells[2],
g->max_nbrs, g->max_atoms,
(int)
((g->total*sizeof(grid_cell)+g->total*g->max_nbrs*sizeof(int*) +
g->total*g->max_nbrs*sizeof(rvec) +
(g->native_end[0]-g->native_str[0])*
(g->native_end[1]-g->native_str[1])*
(g->native_end[2]-g->native_str[2])*g->max_atoms*sizeof(int))/
(1024*1024)) );
#endif
}
@ -691,15 +530,6 @@ void Deallocate_Grid( grid *g )
}
/************ mpi buffers ********************/
/* make the allocations based on the largest need by the three comms:
1- transfer an atom who has moved into other proc's domain (mpi_atom)
2- exchange boundary atoms (boundary_atom)
3- update position info for boundary atoms (mpi_rvec)
the largest space by far is required for the 2nd comm operation above.
buffers are void*, type cast to the correct pointer type to access
the allocated buffers */
int Allocate_MPI_Buffers( mpi_datatypes *mpi_data, int est_recv,
neighbor_proc *my_nbrs, char *msg )
{
@ -766,26 +596,6 @@ void ReAllocate( reax_system *system, control_params *control,
g = &(system->my_grid);
comm = mpi_data->world;
#if defined(DEBUG)
fprintf( stderr, "p%d@reallocate: n: %d, N: %d, numH: %d\n",
system->my_rank, system->n, system->N, system->numH );
fprintf( stderr, "p%d@reallocate: local_cap: %d, total_cap: %d, Hcap: %d\n",
system->my_rank, system->local_cap, system->total_cap,
system->Hcap);
fprintf( stderr, "p%d: realloc.num_far: %d\n",
system->my_rank, realloc->num_far );
fprintf( stderr, "p%d: realloc.H: %d, realloc.Htop: %d\n",
system->my_rank, realloc->H, realloc->Htop );
fprintf( stderr, "p%d: realloc.Hbonds: %d, realloc.num_hbonds: %d\n",
system->my_rank, realloc->hbonds, realloc->num_hbonds );
fprintf( stderr, "p%d: realloc.bonds: %d, num_bonds: %d\n",
system->my_rank, realloc->bonds, realloc->num_bonds );
fprintf( stderr, "p%d: realloc.num_3body: %d\n",
system->my_rank, realloc->num_3body );
#endif
// IMPORTANT: LOOSE ZONES CHECKS ARE DISABLED FOR NOW BY &&'ing with 0!!!
int nflag = 0;
if( system->n >= DANGER_ZONE * system->local_cap ||
(0 && system->n <= LOOSE_ZONE * system->local_cap) ) {
@ -802,12 +612,6 @@ void ReAllocate( reax_system *system, control_params *control,
if( Nflag ) {
/* system */
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating system and workspace -"\
"n=%d N=%d local_cap=%d total_cap=%d\n",
system->my_rank, system->n, system->N,
system->local_cap, system->total_cap );
#endif
ret = Allocate_System( system, system->local_cap, system->total_cap, msg );
if( ret != SUCCESS ) {
fprintf( stderr, "not enough space for atom_list: total_cap=%d",
@ -843,46 +647,12 @@ void ReAllocate( reax_system *system, control_params *control,
newsize = static_cast<int>
(MAX( realloc->num_far*safezone, mincap*MIN_NBRS ));
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating far_nbrs: num_fars=%d, space=%dMB\n",
system->my_rank, (int)(realloc->num_far*SAFE_ZONE),
(newsize*sizeof(far_neighbor_data)/(1024*1024)) );
#endif
Reallocate_Neighbor_List( far_nbrs, system->total_cap, newsize, comm );
realloc->num_far = 0;
}
}
#if defined(PURE_REAX)
/* qeq coef matrix */
H = workspace->H;
if( nflag || realloc->Htop >= H->m * DANGER_ZONE ) {
if( realloc->Htop > H->m ) {
fprintf( stderr,
"step%d - ran out of space on H matrix: Htop=%d, max = %d",
data->step, realloc->Htop, H->m );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating H matrix: Htop=%d, space=%dMB\n",
system->my_rank, (int)(realloc->Htop*SAFE_ZONE),
(int)(realloc->Htop * SAFE_ZONE * sizeof(sparse_matrix_entry) /
(1024*1024)) );
#endif
newsize = static_cast<int>
(MAX( realloc->Htop*safezone, mincap*MIN_NBRS ));
Reallocate_Matrix( &(workspace->H), system->local_cap,
newsize, "H", comm );
//Deallocate_Matrix( workspace->L );
//Deallocate_Matrix( workspace->U );
workspace->L = NULL;
workspace->U = NULL;
realloc->Htop = 0;
}
#endif /*PURE_REAX*/
/* hydrogen bonds list */
if( control->hbond_cut > 0 ) {
Hflag = 0;
@ -895,10 +665,6 @@ void ReAllocate( reax_system *system, control_params *control,
if( Hflag || realloc->hbonds ) {
ret = Reallocate_HBonds_List( system, (*lists)+HBONDS, comm );
realloc->hbonds = 0;
#if defined(DEBUG_FOCUS)
fprintf(stderr, "p%d: reallocating hbonds: total_hbonds=%d space=%dMB\n",
system->my_rank, ret, (int)(ret*sizeof(hbond_data)/(1024*1024)));
#endif
}
}
@ -909,21 +675,10 @@ void ReAllocate( reax_system *system, control_params *control,
&est_3body, comm );
realloc->bonds = 0;
realloc->num_3body = MAX( realloc->num_3body, est_3body );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating bonds: total_bonds=%d, space=%dMB\n",
system->my_rank, num_bonds,
(int)(num_bonds*sizeof(bond_data)/(1024*1024)) );
#endif
}
/* 3-body list */
if( realloc->num_3body > 0 ) {
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating 3body list: num_3body=%d, space=%dMB\n",
system->my_rank, realloc->num_3body,
(int)(realloc->num_3body * sizeof(three_body_interaction_data) /
(1024*1024)) );
#endif
Delete_List( (*lists)+THREE_BODIES, comm );
if( num_bonds == -1 )
@ -939,88 +694,4 @@ void ReAllocate( reax_system *system, control_params *control,
realloc->num_3body = -1;
}
#if defined(PURE_REAX)
/* grid */
if( renbr && realloc->gcell_atoms > -1 ) {
#if defined(DEBUG_FOCUS)
fprintf(stderr, "reallocating gcell: g->max_atoms: %d\n", g->max_atoms);
#endif
for( i = g->native_str[0]; i < g->native_end[0]; i++ )
for( j = g->native_str[1]; j < g->native_end[1]; j++ )
for( k = g->native_str[2]; k < g->native_end[2]; k++ ) {
// reallocate g->atoms
sfree( g->cells[i][j][k].atoms, "g:atoms" );
g->cells[i][j][k].atoms = (int*)
scalloc( realloc->gcell_atoms, sizeof(int), "g:atoms", comm);
}
realloc->gcell_atoms = -1;
}
/* mpi buffers */
// we have to be at a renbring step -
// to ensure correct values at mpi_buffers for update_boundary_positions
if( !renbr )
mpi_flag = 0;
// check whether in_buffer capacity is enough
else if( system->max_recved >= system->est_recv * 0.90 )
mpi_flag = 1;
else {
// otherwise check individual outgoing buffers
mpi_flag = 0;
for( p = 0; p < MAX_NBRS; ++p ) {
nbr_pr = &( system->my_nbrs[p] );
nbr_data = &( mpi_data->out_buffers[p] );
if( nbr_data->cnt >= nbr_pr->est_send * 0.90 ) {
mpi_flag = 1;
break;
}
}
}
if( mpi_flag ) {
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating mpi_buf: old_recv=%d\n",
system->my_rank, system->est_recv );
for( p = 0; p < MAX_NBRS; ++p )
fprintf( stderr, "p%d: nbr%d old_send=%d\n",
system->my_rank, p, system->my_nbrs[p].est_send );
#endif
/* update mpi buffer estimates based on last comm */
system->est_recv = MAX( system->max_recved*SAFER_ZONE, MIN_SEND );
system->est_trans =
(system->est_recv * sizeof(boundary_atom)) / sizeof(mpi_atom);
total_send = 0;
for( p = 0; p < MAX_NBRS; ++p ) {
nbr_pr = &( system->my_nbrs[p] );
nbr_data = &( mpi_data->out_buffers[p] );
nbr_pr->est_send = MAX( nbr_data->cnt*SAFER_ZONE, MIN_SEND );
total_send += nbr_pr->est_send;
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating mpi_buf: recv=%d send=%d total=%dMB\n",
system->my_rank, system->est_recv, total_send,
(int)((system->est_recv+total_send)*sizeof(boundary_atom)/
(1024*1024)));
for( p = 0; p < MAX_NBRS; ++p )
fprintf( stderr, "p%d: nbr%d new_send=%d\n",
system->my_rank, p, system->my_nbrs[p].est_send );
#endif
/* reallocate mpi buffers */
Deallocate_MPI_Buffers( mpi_data );
ret = Allocate_MPI_Buffers( mpi_data, system->est_recv,
system->my_nbrs, msg );
if( ret != SUCCESS ) {
fprintf( stderr, "%s", msg );
fprintf( stderr, "terminating...\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
}
#endif /*PURE_REAX*/
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d @ step%d: reallocate done\n",
system->my_rank, data->step );
MPI_Barrier( comm );
#endif
}

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

@ -25,204 +25,9 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "basic_comm.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_basic_comm.h"
#include "reaxc_vector.h"
#endif
#if defined(PURE_REAX)
void real_packer( void *dummy, mpi_out_data *out_buf )
{
int i;
real *buf = (real*) dummy;
real *out = (real*) out_buf->out_atoms;
for( i = 0; i < out_buf->cnt; ++i )
out[i] = buf[ out_buf->index[i] ];
}
void rvec_packer( void *dummy, mpi_out_data *out_buf )
{
int i;
rvec *buf = (rvec*) dummy;
rvec *out = (rvec*)out_buf->out_atoms;
for( i = 0; i < out_buf->cnt; ++i )
memcpy( out[i], buf[ out_buf->index[i] ], sizeof(rvec) );
}
void rvec2_packer( void *dummy, mpi_out_data *out_buf )
{
int i;
rvec2 *buf = (rvec2*) dummy;
rvec2 *out = (rvec2*) out_buf->out_atoms;
for( i = 0; i < out_buf->cnt; ++i )
memcpy( out[i], buf[ out_buf->index[i] ], sizeof(rvec2) );
}
void Dist( reax_system* system, mpi_datatypes *mpi_data,
void *buf, MPI_Datatype type, int scale, dist_packer pack )
{
int d;
mpi_out_data *out_bufs;
MPI_Comm comm;
MPI_Request req1, req2;
MPI_Status stat1, stat2;
neighbor_proc *nbr1, *nbr2;
#if defined(DEBUG)
fprintf( stderr, "p%d dist: entered\n", system->my_rank );
#endif
comm = mpi_data->comm_mesh3D;
out_bufs = mpi_data->out_buffers;
for( d = 0; d < 3; ++d ) {
/* initiate recvs */
nbr1 = &(system->my_nbrs[2*d]);
if( nbr1->atoms_cnt )
MPI_Irecv( buf + nbr1->atoms_str*scale, nbr1->atoms_cnt, type,
nbr1->rank, 2*d+1,comm, &req1 );
nbr2 = &(system->my_nbrs[2*d+1]);
if( nbr2->atoms_cnt )
MPI_Irecv( buf + nbr2->atoms_str*scale, nbr2->atoms_cnt, type,
nbr2->rank, 2*d, comm, &req2 );
/* send both messages in dimension d */
if( out_bufs[2*d].cnt ) {
pack( buf, out_bufs + (2*d) );
MPI_Send( out_bufs[2*d].out_atoms, out_bufs[2*d].cnt, type,
nbr1->rank, 2*d, comm );
}
if( out_bufs[2*d+1].cnt ) {
pack( buf, out_bufs + (2*d+1) );
MPI_Send( out_bufs[2*d+1].out_atoms, out_bufs[2*d+1].cnt, type,
nbr2->rank, 2*d+1, comm );
}
if( nbr1->atoms_cnt ) MPI_Wait( &req1, &stat1 );
if( nbr2->atoms_cnt ) MPI_Wait( &req2, &stat2 );
}
#if defined(DEBUG)
fprintf( stderr, "p%d dist: done\n", system->my_rank );
#endif
}
void real_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
{
int i;
real *in = (real*) dummy_in;
real *buf = (real*) dummy_buf;
for( i = 0; i < out_buf->cnt; ++i )
buf[ out_buf->index[i] ] += in[i];
}
void rvec_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
{
int i;
rvec *in = (rvec*) dummy_in;
rvec *buf = (rvec*) dummy_buf;
for( i = 0; i < out_buf->cnt; ++i ) {
rvec_Add( buf[ out_buf->index[i] ], in[i] );
#if defined(DEBUG)
fprintf( stderr, "rvec_unpacker: cnt=%d i =%d index[i]=%d\n",
out_buf->cnt, i, out_buf->index[i] );
#endif
}
}
void rvec2_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
{
int i;
rvec2 *in = (rvec2*) dummy_in;
rvec2 *buf = (rvec2*) dummy_buf;
for( i = 0; i < out_buf->cnt; ++i ) {
buf[ out_buf->index[i] ][0] += in[i][0];
buf[ out_buf->index[i] ][1] += in[i][1];
}
}
void Coll( reax_system* system, mpi_datatypes *mpi_data,
void *buf, MPI_Datatype type, int scale, coll_unpacker unpack )
{
int d;
void *in1, *in2;
mpi_out_data *out_bufs;
MPI_Comm comm;
MPI_Request req1, req2;
MPI_Status stat1, stat2;
neighbor_proc *nbr1, *nbr2;
#if defined(DEBUG)
fprintf( stderr, "p%d coll: entered\n", system->my_rank );
#endif
comm = mpi_data->comm_mesh3D;
in1 = mpi_data->in1_buffer;
in2 = mpi_data->in2_buffer;
out_bufs = mpi_data->out_buffers;
for( d = 2; d >= 0; --d ) {
/* initiate recvs */
nbr1 = &(system->my_nbrs[2*d]);
if( out_bufs[2*d].cnt )
MPI_Irecv(in1, out_bufs[2*d].cnt, type, nbr1->rank, 2*d+1, comm, &req1);
nbr2 = &(system->my_nbrs[2*d+1]);
if( out_bufs[2*d+1].cnt )
MPI_Irecv(in2, out_bufs[2*d+1].cnt, type, nbr2->rank, 2*d, comm, &req2);
/* send both messages in dimension d */
if( nbr1->atoms_cnt )
MPI_Send( buf + nbr1->atoms_str*scale, nbr1->atoms_cnt, type,
nbr1->rank, 2*d, comm );
if( nbr2->atoms_cnt )
MPI_Send( buf + nbr2->atoms_str*scale, nbr2->atoms_cnt, type,
nbr2->rank, 2*d+1, comm );
#if defined(DEBUG)
fprintf( stderr, "p%d coll[%d] nbr1: str=%d cnt=%d recv=%d\n",
system->my_rank, d, nbr1->atoms_str, nbr1->atoms_cnt,
out_bufs[2*d].cnt );
fprintf( stderr, "p%d coll[%d] nbr2: str=%d cnt=%d recv=%d\n",
system->my_rank, d, nbr2->atoms_str, nbr2->atoms_cnt,
out_bufs[2*d+1].cnt );
#endif
if( out_bufs[2*d].cnt ) {
MPI_Wait( &req1, &stat1 );
unpack( in1, buf, out_bufs + (2*d) );
}
if( out_bufs[2*d+1].cnt ) {
MPI_Wait( &req2, &stat2 );
unpack( in2, buf, out_bufs + (2*d+1) );
}
}
#if defined(DEBUG)
fprintf( stderr, "p%d coll: done\n", system->my_rank );
#endif
}
#endif /*PURE_REAX*/
/*****************************************************************************/
real Parallel_Norm( real *v, int n, MPI_Comm comm )
{
int i;
@ -237,8 +42,6 @@ real Parallel_Norm( real *v, int n, MPI_Comm comm )
return sqrt( norm_sqr );
}
real Parallel_Dot( real *v1, real *v2, int n, MPI_Comm comm )
{
int i;
@ -253,8 +56,6 @@ real Parallel_Dot( real *v1, real *v2, int n, MPI_Comm comm )
return res;
}
real Parallel_Vector_Acc( real *v, int n, MPI_Comm comm )
{
int i;
@ -268,50 +69,3 @@ real Parallel_Vector_Acc( real *v, int n, MPI_Comm comm )
return res;
}
/*****************************************************************************/
#if defined(TEST_FORCES)
void Coll_ids_at_Master( reax_system *system, storage *workspace,
mpi_datatypes *mpi_data )
{
int i;
rc_tagint *id_list;
MPI_Gather( &system->n, 1, MPI_INT, workspace->rcounts, 1, MPI_INT,
MASTER_NODE, mpi_data->world );
if( system->my_rank == MASTER_NODE ){
workspace->displs[0] = 0;
for( i = 1; i < system->wsize; ++i )
workspace->displs[i] = workspace->displs[i-1] + workspace->rcounts[i-1];
}
id_list = (int*) malloc( system->n * sizeof(int) );
for( i = 0; i < system->n; ++i )
id_list[i] = system->my_atoms[i].orig_id;
MPI_Gatherv( id_list, system->n, MPI_INT,
workspace->id_all, workspace->rcounts, workspace->displs,
MPI_INT, MASTER_NODE, mpi_data->world );
free( id_list );
#if defined(DEBUG)
if( system->my_rank == MASTER_NODE ) {
for( i = 0 ; i < system->bigN; ++i )
fprintf( stderr, "id_all[%d]: %d\n", i, workspace->id_all[i] );
}
#endif
}
void Coll_rvecs_at_Master( reax_system *system, storage *workspace,
mpi_datatypes *mpi_data, rvec* v )
{
MPI_Gatherv( v, system->n, mpi_data->mpi_rvec,
workspace->f_all, workspace->rcounts, workspace->displs,
mpi_data->mpi_rvec, MASTER_NODE, mpi_data->world );
}
#endif

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

@ -26,370 +26,9 @@
#include "pair_reax_c.h"
#include "reaxc_types.h"
#if defined(PURE_REAX)
#include "bond_orders.h"
#include "list.h"
#include "vector.h"
#include "io_tools.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_bond_orders.h"
#include "reaxc_list.h"
#include "reaxc_vector.h"
#endif
#ifdef TEST_FORCES
void Get_dBO( reax_system *system, reax_list **lists,
int i, int pj, real C, rvec *v )
{
reax_list *bonds = (*lists) + BONDS;
reax_list *dBOs = (*lists) + DBOS;
int start_pj, end_pj, k;
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
for( k = start_pj; k < end_pj; ++k )
rvec_Scale( v[dBOs->select.dbo_list[k].wrt],
C, dBOs->select.dbo_list[k].dBO );
}
void Get_dBOpinpi2( reax_system *system, reax_list **lists,
int i, int pj, real Cpi, real Cpi2, rvec *vpi, rvec *vpi2 )
{
reax_list *bonds = (*lists) + BONDS;
reax_list *dBOs = (*lists) + DBOS;
dbond_data *dbo_k;
int start_pj, end_pj, k;
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
for( k = start_pj; k < end_pj; ++k )
{
dbo_k = &(dBOs->select.dbo_list[k]);
rvec_Scale( vpi[dbo_k->wrt], Cpi, dbo_k->dBOpi );
rvec_Scale( vpi2[dbo_k->wrt], Cpi2, dbo_k->dBOpi2 );
}
}
void Add_dBO( reax_system *system, reax_list **lists,
int i, int pj, real C, rvec *v )
{
reax_list *bonds = (*lists) + BONDS;
reax_list *dBOs = (*lists) + DBOS;
int start_pj, end_pj, k;
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
//fprintf( stderr, "i=%d j=%d start=%d end=%d\n", i, pj, start_pj, end_pj );
for( k = start_pj; k < end_pj; ++k )
rvec_ScaledAdd( v[dBOs->select.dbo_list[k].wrt],
C, dBOs->select.dbo_list[k].dBO );
}
void Add_dBOpinpi2( reax_system *system, reax_list **lists,
int i, int pj, real Cpi, real Cpi2, rvec *vpi, rvec *vpi2 )
{
reax_list *bonds = (*lists) + BONDS;
reax_list *dBOs = (*lists) + DBOS;
dbond_data *dbo_k;
int start_pj, end_pj, k;
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
for( k = start_pj; k < end_pj; ++k )
{
dbo_k = &(dBOs->select.dbo_list[k]);
rvec_ScaledAdd( vpi[dbo_k->wrt], Cpi, dbo_k->dBOpi );
rvec_ScaledAdd( vpi2[dbo_k->wrt], Cpi2, dbo_k->dBOpi2 );
}
}
void Add_dBO_to_Forces( reax_system *system, reax_list **lists,
int i, int pj, real C )
{
reax_list *bonds = (*lists) + BONDS;
reax_list *dBOs = (*lists) + DBOS;
int start_pj, end_pj, k;
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
for( k = start_pj; k < end_pj; ++k )
rvec_ScaledAdd( system->my_atoms[dBOs->select.dbo_list[k].wrt].f,
C, dBOs->select.dbo_list[k].dBO );
}
void Add_dBOpinpi2_to_Forces( reax_system *system, reax_list **lists,
int i, int pj, real Cpi, real Cpi2 )
{
reax_list *bonds = (*lists) + BONDS;
reax_list *dBOs = (*lists) + DBOS;
dbond_data *dbo_k;
int start_pj, end_pj, k;
pj = bonds->select.bond_list[pj].dbond_index;
start_pj = Start_Index(pj, dBOs);
end_pj = End_Index(pj, dBOs);
for( k = start_pj; k < end_pj; ++k )
{
dbo_k = &(dBOs->select.dbo_list[k]);
rvec_ScaledAdd( system->my_atoms[dbo_k->wrt].f, Cpi, dbo_k->dBOpi );
rvec_ScaledAdd( system->my_atoms[dbo_k->wrt].f, Cpi2, dbo_k->dBOpi2 );
}
}
void Add_dDelta( reax_system *system, reax_list **lists, int i, real C, rvec *v )
{
reax_list *dDeltas = &((*lists)[DDELTAS]);
int start = Start_Index(i, dDeltas);
int end = End_Index(i, dDeltas);
int k;
for( k = start; k < end; ++k )
rvec_ScaledAdd( v[dDeltas->select.dDelta_list[k].wrt],
C, dDeltas->select.dDelta_list[k].dVal );
}
void Add_dDelta_to_Forces( reax_system *system, reax_list **lists,
int i, real C )
{
reax_list *dDeltas = &((*lists)[DDELTAS]);
int start = Start_Index(i, dDeltas);
int end = End_Index(i, dDeltas);
int k;
for( k = start; k < end; ++k )
rvec_ScaledAdd( system->my_atoms[dDeltas->select.dDelta_list[k].wrt].f,
C, dDeltas->select.dDelta_list[k].dVal );
}
void Calculate_dBO( int i, int pj,
storage *workspace, reax_list **lists, int *top )
{
/* Initializations */
reax_list *bonds, *dBOs;
int j, k, l, start_i, end_i, end_j;
bond_data *nbr_l, *nbr_k;
bond_order_data *bo_ij;
dbond_data *top_dbo;
// rvec due_j[1000], due_i[1000];
// rvec due_j_pi[1000], due_i_pi[1000];
// memset(due_j, 0, sizeof(rvec)*1000 );
// memset(due_i, 0, sizeof(rvec)*1000 );
// memset(due_j_pi, 0, sizeof(rvec)*1000 );
// memset(due_i_pi, 0, sizeof(rvec)*1000 );
bonds = *lists + BONDS;
dBOs = *lists + DBOS;
start_i = Start_Index(i, bonds);
end_i = End_Index(i, bonds);
j = bonds->select.bond_list[pj].nbr;
l = Start_Index(j, bonds);
end_j = End_Index(j, bonds);
bo_ij = &(bonds->select.bond_list[pj].bo_data);
top_dbo = &(dBOs->select.dbo_list[ (*top) ]);
for( k = start_i; k < end_i; ++k ) {
nbr_k = &(bonds->select.bond_list[k]);
for( ; l < end_j && bonds->select.bond_list[l].nbr < nbr_k->nbr; ++l ) {
/* These are the neighbors of j which are not in the nbr_list of i
Note that they might also include i! */
nbr_l = &(bonds->select.bond_list[l]);
top_dbo->wrt = nbr_l->nbr;
//rvec_ScaledAdd( due_j[top_dbo->wrt],
// -bo_ij->BO * bo_ij->A2_ji, nbr_l->bo_data.dBOp );
/*3rd, dBO*/
rvec_Scale( top_dbo->dBO, -bo_ij->C3dbo, nbr_l->bo_data.dBOp );
/*4th, dBOpi*/
rvec_Scale( top_dbo->dBOpi, -bo_ij->C4dbopi, nbr_l->bo_data.dBOp );
/*4th, dBOpp*/
rvec_Scale( top_dbo->dBOpi2, -bo_ij->C4dbopi2, nbr_l->bo_data.dBOp );
if( nbr_l->nbr == i ) {
/* do the adjustments on i */
//rvec_ScaledAdd( due_i[i], bo_ij->A0_ij+bo_ij->BO*bo_ij->A1_ij,
// bo_ij->dBOp ); /*1st, dBO*/
//rvec_ScaledAdd( due_i[i], bo_ij->BO * bo_ij->A2_ij,
// workspace->dDeltap_self[i] ); /*2nd, dBO*/
/*1st, dBO*/
rvec_ScaledAdd( top_dbo->dBO, bo_ij->C1dbo, bo_ij->dBOp );
/*2nd, dBO*/
rvec_ScaledAdd( top_dbo->dBO,
bo_ij->C2dbo, workspace->dDeltap_self[i] );
/*1st, dBOpi*/
rvec_ScaledAdd(top_dbo->dBOpi,bo_ij->C1dbopi,bo_ij->dln_BOp_pi );
/*2nd, dBOpi*/
rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C2dbopi, bo_ij->dBOp );
/*3rd, dBOpi*/
rvec_ScaledAdd( top_dbo->dBOpi,
bo_ij->C3dbopi, workspace->dDeltap_self[i] );
/*1st, dBO_p*/
rvec_ScaledAdd( top_dbo->dBOpi2,
bo_ij->C1dbopi2, bo_ij->dln_BOp_pi2 );
/*2nd, dBO_p*/
rvec_ScaledAdd( top_dbo->dBOpi2, bo_ij->C2dbopi2, bo_ij->dBOp );
/*3rd, dBO_p*/
rvec_ScaledAdd( top_dbo->dBOpi2,
bo_ij->C3dbopi2, workspace->dDeltap_self[i] );
}
//rvec_Add( workspace->dDelta[nbr_l->nbr], top_dbo->dBO );
++(*top), ++top_dbo;
}
/* Now we are processing neighbor k of i. */
top_dbo->wrt = nbr_k->nbr;
//2nd, dBO
//rvec_ScaledAdd( due_i[top_dbo->wrt],
// -bo_ij->BO * bo_ij->A2_ij, nbr_k->bo_data.dBOp );
/*2nd, dBO*/
rvec_Scale( top_dbo->dBO, -bo_ij->C2dbo, nbr_k->bo_data.dBOp );
/*3rd, dBOpi*/
rvec_Scale( top_dbo->dBOpi, -bo_ij->C3dbopi, nbr_k->bo_data.dBOp );
/*3rd, dBOpp*/
rvec_Scale( top_dbo->dBOpi2, -bo_ij->C3dbopi2, nbr_k->bo_data.dBOp );
if( l < end_j && bonds->select.bond_list[l].nbr == nbr_k->nbr ) {
/* This is a common neighbor of i and j. */
nbr_l = &(bonds->select.bond_list[l]);
/*3rd, dBO*/
//rvec_ScaledAdd( due_j[top_dbo->wrt],
// -bo_ij->BO * bo_ij->A2_ji, nbr_l->bo_data.dBOp );
/*3rd, dBO*/
rvec_ScaledAdd( top_dbo->dBO, -bo_ij->C3dbo, nbr_l->bo_data.dBOp );
/*4th, dBOpi*/
rvec_ScaledAdd(top_dbo->dBOpi, -bo_ij->C4dbopi, nbr_l->bo_data.dBOp);
/*4th, dBOpp*/
rvec_ScaledAdd(top_dbo->dBOpi2,-bo_ij->C4dbopi2,nbr_l->bo_data.dBOp);
++l;
}
else if( k == pj ) {
/* This negihbor is j. */
//rvec_ScaledAdd( due_j[j], -(bo_ij->A0_ij+bo_ij->BO*bo_ij->A1_ij),
// bo_ij->dBOp ); /*1st, dBO*/
//rvec_ScaledAdd( due_j[j], bo_ij->BO * bo_ij->A2_ji,
// workspace->dDeltap_self[j] ); /*3rd, dBO*/
/*1st, dBO*/
rvec_ScaledAdd( top_dbo->dBO, -bo_ij->C1dbo, bo_ij->dBOp );
/*3rd, dBO*/
rvec_ScaledAdd(top_dbo->dBO,bo_ij->C3dbo,workspace->dDeltap_self[j]);
/*1st, dBOpi*/
rvec_ScaledAdd( top_dbo->dBOpi, -bo_ij->C1dbopi, bo_ij->dln_BOp_pi );
/*2nd, dBOpi*/
rvec_ScaledAdd( top_dbo->dBOpi, -bo_ij->C2dbopi, bo_ij->dBOp );
/*4th, dBOpi*/
rvec_ScaledAdd( top_dbo->dBOpi,
bo_ij->C4dbopi, workspace->dDeltap_self[j] );
/*1st, dBOpi2*/
rvec_ScaledAdd(top_dbo->dBOpi2,-bo_ij->C1dbopi2,bo_ij->dln_BOp_pi2);
/*2nd, dBOpi2*/
rvec_ScaledAdd(top_dbo->dBOpi2,-bo_ij->C2dbopi2,bo_ij->dBOp );
/*4th, dBOpi2*/
rvec_ScaledAdd(top_dbo->dBOpi2, bo_ij->C4dbopi2,
workspace->dDeltap_self[j] );
}
// rvec_Add( workspace->dDelta[nbr_k->nbr], top_dbo->dBO );
++(*top), ++top_dbo;
}
for( ; l < end_j; ++l ) {
/* These are the remaining neighbors of j which are not in the
neighbor_list of i. Note that they might also include i!*/
nbr_l = &(bonds->select.bond_list[l]);
top_dbo->wrt = nbr_l->nbr;
// fprintf( stdout, "\tl: %d nbr:%d\n", l, nbr_l->nbr+1 );
// rvec_ScaledAdd( due_j[top_dbo->wrt],
// -bo_ij->BO * bo_ij->A2_ji, nbr_l->bo_data.dBOp );
/*3rd, dBO*/
rvec_Scale( top_dbo->dBO, -bo_ij->C3dbo, nbr_l->bo_data.dBOp );
/*4th, dBOpi*/
rvec_Scale( top_dbo->dBOpi, -bo_ij->C4dbopi, nbr_l->bo_data.dBOp );
/*4th, dBOpp*/
rvec_Scale( top_dbo->dBOpi2, -bo_ij->C4dbopi2, nbr_l->bo_data.dBOp );
if( nbr_l->nbr == i )
{
/* do the adjustments on i */
//rvec_ScaledAdd( due_i[i], bo_ij->A0_ij + bo_ij->BO * bo_ij->A1_ij,
// bo_ij->dBOp ); /*1st, dBO*/
//rvec_ScaledAdd( due_i[i], bo_ij->BO * bo_ij->A2_ij,
// workspace->dDeltap_self[i] ); /*2nd, dBO*/
/*1st, dBO*/
rvec_ScaledAdd( top_dbo->dBO, bo_ij->C1dbo, bo_ij->dBOp );
/*2nd, dBO*/
rvec_ScaledAdd(top_dbo->dBO,bo_ij->C2dbo,workspace->dDeltap_self[i]);
/*1st, dBO_p*/
rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C1dbopi, bo_ij->dln_BOp_pi );
/*2nd, dBOpi*/
rvec_ScaledAdd( top_dbo->dBOpi, bo_ij->C2dbopi, bo_ij->dBOp );
/*3rd,dBOpi*/
rvec_ScaledAdd( top_dbo->dBOpi,
bo_ij->C3dbopi, workspace->dDeltap_self[i] );
/*1st, dBO_p*/
rvec_ScaledAdd(top_dbo->dBOpi2, bo_ij->C1dbopi2, bo_ij->dln_BOp_pi2);
/*2nd, dBO_p*/
rvec_ScaledAdd(top_dbo->dBOpi2, bo_ij->C2dbopi2, bo_ij->dBOp);
/*3rd,dBO_p*/
rvec_ScaledAdd(top_dbo->dBOpi2,
bo_ij->C3dbopi2, workspace->dDeltap_self[i]);
}
// rvec_Add( workspace->dDelta[nbr_l->nbr], top_dbo->dBO );
++(*top), ++top_dbo;
}
// for( k = 0; k < 21; ++k ){
// fprintf( stderr, "%d %d %d, due_i:[%g %g %g]\n",
// i+1, j+1, k+1, due_i[k][0], due_i[k][1], due_i[k][2] );
// fprintf( stderr, "%d %d %d, due_j:[%g %g %g]\n",
// i+1, j+1, k+1, due_j[k][0], due_j[k][1], due_j[k][2] );
}
#endif
void Add_dBond_to_Forces_NPT( int i, int pj, simulation_data *data,
storage *workspace, reax_list **lists )
@ -426,11 +65,6 @@ void Add_dBond_to_Forces_NPT( int i, int pj, simulation_data *data,
coef.C2dDelta = bo_ij->C2dbo * (workspace->CdDelta[i]+workspace->CdDelta[j]);
coef.C3dDelta = bo_ij->C3dbo * (workspace->CdDelta[i]+workspace->CdDelta[j]);
/************************************
* forces related to atom i *
* first neighbors of atom i *
************************************/
for( pk = Start_Index(i, bonds); pk < End_Index(i, bonds); ++pk ) {
nbr_k = &(bonds->select.bond_list[pk]);
k = nbr_k->nbr;
@ -446,16 +80,6 @@ void Add_dBond_to_Forces_NPT( int i, int pj, simulation_data *data,
rvec_iMultiply( ext_press, nbr_k->rel_box, temp );
rvec_Add( data->my_ext_press, ext_press );
/* if( !ivec_isZero( nbr_k->rel_box ) )
fprintf( stderr, "%3d %3d %3d: dvec[%10.6f %10.6f %10.6f]"
"ext[%3d %3d %3d] f[%10.6f %10.6f %10.6f]\n",
i+1, system->my_atoms[i].x[0],
system->my_atoms[i].x[1], system->my_atoms[i].x[2],
j+1, k+1, system->my_atoms[k].x[0],
system->my_atoms[k].x[1], system->my_atoms[k].x[2],
nbr_k->dvec[0], nbr_k->dvec[1], nbr_k->dvec[2],
nbr_k->rel_box[0], nbr_k->rel_box[1], nbr_k->rel_box[2],
temp[0], temp[1], temp[2] ); */
}
/* then atom i itself */
@ -473,13 +97,7 @@ void Add_dBond_to_Forces_NPT( int i, int pj, simulation_data *data,
/* force */
rvec_Add( workspace->f[i], temp );
/* ext pressure due to i is dropped, counting force on j will be enough */
/******************************************************
* forces and pressure related to atom j *
* first neighbors of atom j *
******************************************************/
for( pk = Start_Index(j, bonds); pk < End_Index(j, bonds); ++pk ) {
nbr_k = &(bonds->select.bond_list[pk]);
k = nbr_k->nbr;
@ -497,16 +115,6 @@ void Add_dBond_to_Forces_NPT( int i, int pj, simulation_data *data,
rvec_iMultiply( ext_press, rel_box, temp );
rvec_Add( data->my_ext_press, ext_press );
/* if( !ivec_isZero( rel_box ) )
fprintf( stderr, "%3d %3d %3d: dvec[%10.6f %10.6f %10.6f]"
"ext[%3d %3d %3d] f[%10.6f %10.6f %10.6f]\n",
i+1, j+1, system->my_atoms[j].x[0],
system->my_atoms[j].x[1], system->my_atoms[j].x[2],
k+1, system->my_atoms[k].x[0],
system->my_atoms[k].x[1], system->my_atoms[k].x[2],
nbr_k->dvec[0], nbr_k->dvec[1], nbr_k->dvec[2],
rel_box[0], rel_box[1], rel_box[2],
temp[0], temp[1], temp[2] ); */
}
}
@ -530,20 +138,8 @@ void Add_dBond_to_Forces_NPT( int i, int pj, simulation_data *data,
rvec_iMultiply( ext_press, nbr_j->rel_box, temp );
rvec_Add( data->my_ext_press, ext_press );
/* if( !ivec_isZero( nbr_j->rel_box ) )
fprintf( stderr, "%3d %3d %3d: dvec[%10.6f %10.6f %10.6f]"
"ext[%3d %3d %3d] f[%10.6f %10.6f %10.6f]\n",
i+1, system->my_atoms[i].x[0], system->my_atoms[i].x[1],
system->my_atoms[i].x[2],
j+1,system->my_atoms[j].x[0], system->my_atoms[j].x[1],
system->my_atoms[j].x[2],
j+1, nbr_j->dvec[0], nbr_j->dvec[1], nbr_j->dvec[2],
nbr_j->rel_box[0], nbr_j->rel_box[1], nbr_j->rel_box[2],
temp[0], temp[1], temp[2] ); */
}
void Add_dBond_to_Forces( reax_system *system, int i, int pj,
storage *workspace, reax_list **lists )
{
@ -554,6 +150,7 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
int pk, k, j;
/* Virial Tallying variables */
real f_scaler;
rvec fi_tmp, fj_tmp, fk_tmp, delij, delji, delki, delkj, temp;
/* Initializations */
@ -738,8 +335,6 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
rvec_Scale(bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi );
rvec_Scale(bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2 );
/* Only dBOp wrt. dr_i is stored here, note that
dBOp/dr_i = -dBOp/dr_j and all others are 0 */
rvec_Scale( bo_ij->dBOp,
-(bo_ij->BO_s * Cln_BOp_s +
bo_ij->BO_pi * Cln_BOp_pi +
@ -758,35 +353,6 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
bo_ij->Cdbo = bo_ij->Cdbopi = bo_ij->Cdbopi2 = 0.0;
bo_ji->Cdbo = bo_ji->Cdbopi = bo_ji->Cdbopi2 = 0.0;
/*fprintf( stderr, "%d %d %g %g %g\n",
i+1, j+1, bo_ij->BO, bo_ij->BO_pi, bo_ij->BO_pi2 );*/
/*fprintf( stderr, "Cln_BOp_s: %f, pbo2: %f, C12:%f\n",
Cln_BOp_s, twbp->p_bo2, C12 );
fprintf( stderr, "Cln_BOp_pi: %f, pbo4: %f, C34:%f\n",
Cln_BOp_pi, twbp->p_bo4, C34 );
fprintf( stderr, "Cln_BOp_pi2: %f, pbo6: %f, C56:%f\n",
Cln_BOp_pi2, twbp->p_bo6, C56 );*/
/*fprintf(stderr, "pbo1: %f, pbo2:%f\n", twbp->p_bo1, twbp->p_bo2);
fprintf(stderr, "pbo3: %f, pbo4:%f\n", twbp->p_bo3, twbp->p_bo4);
fprintf(stderr, "pbo5: %f, pbo6:%f\n", twbp->p_bo5, twbp->p_bo6);
fprintf( stderr, "r_s: %f, r_p: %f, r_pp: %f\n",
twbp->r_s, twbp->r_p, twbp->r_pp );
fprintf( stderr, "C12: %g, C34:%g, C56:%g\n", C12, C34, C56 );*/
/*fprintf( stderr, "\tfactors: %g %g %g\n",
-(bo_ij->BO_s * Cln_BOp_s + bo_ij->BO_pi * Cln_BOp_pi +
bo_ij->BO_pi2 * Cln_BOp_pp),
-bo_ij->BO_pi * Cln_BOp_pi, -bo_ij->BO_pi2 * Cln_BOp_pi2 );*/
/*fprintf( stderr, "dBOpi:\t[%g, %g, %g]\n",
bo_ij->dBOp[0], bo_ij->dBOp[1], bo_ij->dBOp[2] );
fprintf( stderr, "dBOpi:\t[%g, %g, %g]\n",
bo_ij->dln_BOp_pi[0], bo_ij->dln_BOp_pi[1],
bo_ij->dln_BOp_pi[2] );
fprintf( stderr, "dBOpi2:\t[%g, %g, %g]\n\n",
bo_ij->dln_BOp_pi2[0], bo_ij->dln_BOp_pi2[1],
bo_ij->dln_BOp_pi2[2] );*/
return 1;
}
@ -804,7 +370,7 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
storage *workspace, reax_list **lists, output_controls *out_control )
{
int i, j, pj, type_i, type_j;
int start_i, end_i, sym_index;
int start_i, end_i, sym_index, num_bonds;
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;
@ -817,71 +383,45 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
two_body_parameters *twbp;
bond_order_data *bo_ij, *bo_ji;
reax_list *bonds = (*lists) + BONDS;
#ifdef TEST_FORCES
int k, pk, start_j, end_j;
int top_dbo, top_dDelta;
dbond_data *pdbo;
dDelta_data *ptop_dDelta;
reax_list *dDeltas, *dBOs;
top_dbo=0;
top_dDelta=0;
dDeltas = (*lists) + DDELTAS;
dBOs = (*lists) + DBOS;
//for( i = 0; i < system->N; ++i )
// qsort( &(bonds->select.bond_list[Start_Index(i, bonds)]),
// Num_Entries(i, bonds), sizeof(bond_data), compare_bonds );
#endif
num_bonds = 0;
p_boc1 = system->reax_param.gp.l[0];
p_boc2 = system->reax_param.gp.l[1];
/* Calculate Deltaprime, Deltaprime_boc values */
for( i = 0; i < system->N; ++i ) {
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
sbp_i = &(system->reax_param.sbp[type_i]);
workspace->Deltap[i] = workspace->total_bond_order[i] - sbp_i->valency;
workspace->Deltap_boc[i] =
workspace->total_bond_order[i] - sbp_i->valency_val;
//fprintf( stdout, "%d(%d) %24.15f\n",
// i, workspace->bond_mark[i], workspace->total_bond_order[i] );
workspace->total_bond_order[i] = 0;
}
// fprintf( stderr, "done with uncorrected bond orders\n" );
/* Corrected Bond Order calculations */
for( i = 0; i < system->N; ++i ) {
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
sbp_i = &(system->reax_param.sbp[type_i]);
val_i = sbp_i->valency;
Deltap_i = workspace->Deltap[i];
Deltap_boc_i = workspace->Deltap_boc[i];
start_i = Start_Index(i, bonds);
end_i = End_Index(i, bonds);
// fprintf( stderr, "i:%d Dp:%g Dbocp:%g s:%d e:%d\n",
// i+1, Deltap_i, Deltap_boc_i, start_i, end_i );
for( pj = start_i; pj < end_i; ++pj ) {
j = bonds->select.bond_list[pj].nbr;
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
bo_ij = &( bonds->select.bond_list[pj].bo_data );
// fprintf( stderr, "\tj:%d - ubo: %8.3f\n", j+1, bo_ij->BO );
if( i < j || workspace->bond_mark[j] > 3 ) {
twbp = &( system->reax_param.tbp[type_i][type_j] );
#ifdef TEST_FORCES
Set_Start_Index( pj, top_dbo, dBOs );
/* fprintf( stderr, "%6d%6d%12.6f%12.6f%12.6f\n",
workspace->reverse_map[i], workspace->reverse_map[j],
twbp->ovc, twbp->v13cor, bo_ij->BO ); */
#endif
if( twbp->ovc < 0.001 && twbp->v13cor < 0.001 ) {
/* There is no correction to bond orders nor to derivatives
of bond order prime! So we leave bond orders unchanged and
set derivative of bond order coefficients such that
dBO = dBOp & dBOxx = dBOxxp in Add_dBO_to_Forces */
bo_ij->C1dbo = 1.000000;
bo_ij->C2dbo = 0.000000;
bo_ij->C3dbo = 0.000000;
@ -896,24 +436,6 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
bo_ij->C3dbopi2 = 0.000000;
bo_ij->C4dbopi2 = 0.000000;
#ifdef TEST_FORCES
pdbo = &(dBOs->select.dbo_list[ top_dbo ]);
// compute dBO_ij/dr_i
pdbo->wrt = i;
rvec_Copy( pdbo->dBO, bo_ij->dBOp );
rvec_Scale( pdbo->dBOpi, bo_ij->BO_pi, bo_ij->dln_BOp_pi );
rvec_Scale( pdbo->dBOpi2, bo_ij->BO_pi2, bo_ij->dln_BOp_pi2);
// compute dBO_ij/dr_j
pdbo++;
pdbo->wrt = j;
rvec_Scale( pdbo->dBO, -1.0, bo_ij->dBOp );
rvec_Scale( pdbo->dBOpi, -bo_ij->BO_pi, bo_ij->dln_BOp_pi );
rvec_Scale(pdbo->dBOpi2, -bo_ij->BO_pi2, bo_ij->dln_BOp_pi2);
top_dbo += 2;
#endif
}
else {
val_j = system->reax_param.sbp[type_j].valency;
@ -933,14 +455,6 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
f1 = 0.5 * ( ( val_i + f2 )/( val_i + f2 + f3 ) +
( val_j + f2 )/( val_j + f2 + f3 ) );
/*fprintf( stderr,"%d %d\t%g %g j:%g %g p_boc:%g %g\n"
"\tf:%g %g %g, exp:%g %g %g %g\n",
i+1, j+1,
val_i, Deltap_i, val_j, Deltap_j, p_boc1, p_boc2,
f1, f2, f3, exp_p1i, exp_p2i, exp_p1j, exp_p2j );*/
/* Now come the derivates */
/* Bond Order pages 5-7, derivative of f1 */
temp = f2 + f3;
u1_ij = val_i + temp;
u1_ji = val_j + temp;
@ -949,8 +463,6 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
Cf1B_ij = -0.5 * (( u1_ij - f3 ) / SQR( u1_ij ) +
( u1_ji - f3 ) / SQR( u1_ji ));
//Cf1_ij = -Cf1A_ij * p_boc1 * exp_p1i +
// Cf1B_ij * exp_p2i / ( exp_p2i + exp_p2j );
Cf1_ij = 0.50 * ( -p_boc1 * exp_p1i / u1_ij -
((val_i+f2) / SQR(u1_ij)) *
( -p_boc1 * exp_p1i +
@ -964,7 +476,6 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
Cf1_ji = -Cf1A_ij * p_boc1 * exp_p1j +
Cf1B_ij * exp_p2j / ( exp_p2i + exp_p2j );
//fprintf( stderr, "\tCf1:%g %g\n", Cf1_ij, Cf1_ji );
}
else {
/* No overcoordination correction! */
@ -984,12 +495,6 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
f4f5 = f4 * f5;
/* Bond Order pages 8-9, derivative of f4 and f5 */
/*temp = twbp->p_boc5 -
twbp->p_boc3 * twbp->p_boc4 * SQR( bo_ij->BO );
u_ij = temp + twbp->p_boc3 * Deltap_boc_i;
u_ji = temp + twbp->p_boc3 * Deltap_boc_j;
Cf45_ij = Cf45( u_ij, u_ji ) / f4f5;
Cf45_ji = Cf45( u_ji, u_ij ) / f4f5;*/
Cf45_ij = -f4 * exp_f4;
Cf45_ji = -f5 * exp_f5;
}
@ -1007,12 +512,6 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
A3_ij = A2_ij + Cf1_ij / f1;
A3_ji = A2_ji + Cf1_ji / f1;
/*fprintf( stderr, "\tBO: %f, A0: %f, A1: %f"
"A2_ij: %f A2_ji: %f, A3_ij: %f, A3_ji: %f\n",
bo_ij->BO,
A0_ij, A1_ij, A2_ij, A2_ji, A3_ij, A3_ji );*/
/* find corrected bond orders and their derivative coef */
bo_ij->BO = bo_ij->BO * A0_ij;
bo_ij->BO_pi = bo_ij->BO_pi * A0_ij *f1;
@ -1033,9 +532,6 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
bo_ij->C3dbopi2 = bo_ij->BO_pi2 * A3_ij;
bo_ij->C4dbopi2 = bo_ij->BO_pi2 * A3_ji;
#ifdef TEST_FORCES
Calculate_dBO( i, pj, workspace, lists, &top_dbo );
#endif
}
/* neglect bonds that are < 1e-10 */
@ -1050,27 +546,6 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
workspace->total_bond_order[i] += bo_ij->BO; //now keeps total_BO
/* fprintf( stderr, "%d %d\t%g %g %g %g\n"
"Cdbo:\t%g %g %g\n"
"Cdbopi:\t%g %g %g %g\n"
"Cdbopi2:%g %g %g %g\n\n",
i+1, j+1,
bonds->select.bond_list[ pj ].d,
bo_ij->BO,bo_ij->BO_pi, bo_ij->BO_pi2,
bo_ij->C1dbo, bo_ij->C2dbo, bo_ij->C3dbo,
bo_ij->C1dbopi, bo_ij->C2dbopi,
bo_ij->C3dbopi, bo_ij->C4dbopi,
bo_ij->C1dbopi2,bo_ij->C2dbopi2,
bo_ij->C3dbopi2, bo_ij->C4dbopi2 ); */
/* fprintf( stderr, "%d %d BO:%f BO_s:%f BO_pi:%f BO_pi2:%f\n",
i+1,j+1,bo_ij->BO,bo_ij->BO_s,bo_ij->BO_pi,bo_ij->BO_pi2 );*/
#ifdef TEST_FORCES
Set_End_Index( pj, top_dbo, dBOs );
Add_dBO( system, lists, i, pj, 1.0, workspace->dDelta );
#endif
}
else {
/* We only need to update bond orders from bo_ji
@ -1083,63 +558,15 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
bo_ij->BO_pi2 = bo_ji->BO_pi2;
workspace->total_bond_order[i] += bo_ij->BO;// now keeps total_BO
#ifdef TEST_FORCES
Add_dBO( system, lists, j, sym_index, 1.0, workspace->dDelta );
#endif
}
}
#ifdef TEST_FORCES
// fprintf( stderr, "dDelta computations\nj:" );
Set_Start_Index( i, top_dDelta, dDeltas );
ptop_dDelta = &( dDeltas->select.dDelta_list[top_dDelta] );
for( pj = start_i; pj < end_i; ++pj ) {
j = bonds->select.bond_list[pj].nbr;
// fprintf( stderr, "%d ", j );
if( !rvec_isZero( workspace->dDelta[j] ) ) {
ptop_dDelta->wrt = j;
rvec_Copy( ptop_dDelta->dVal, workspace->dDelta[j] );
rvec_MakeZero( workspace->dDelta[j] );
++top_dDelta, ++ptop_dDelta;
}
start_j = Start_Index(j, bonds);
end_j = End_Index(j, bonds);
for( pk = start_j; pk < end_j; ++pk ) {
k = bonds->select.bond_list[pk].nbr;
if( !rvec_isZero( workspace->dDelta[k] ) ) {
ptop_dDelta->wrt = k;
rvec_Copy( ptop_dDelta->dVal, workspace->dDelta[k] );
rvec_MakeZero( workspace->dDelta[k] );
++top_dDelta, ++ptop_dDelta;
}
}
}
Set_End_Index( i, top_dDelta, dDeltas );
/*for(pj = Start_Index(i,dDeltas); pj < End_Index(i,dDeltas); ++pj)
fprintf( stdout, "dDel: %d %d [%g %g %g]\n",
i+1, dDeltas->select.dDelta_list[pj].wrt+1,
dDeltas->select.dDelta_list[pj].dVal[0],
dDeltas->select.dDelta_list[pj].dVal[1],
dDeltas->select.dDelta_list[pj].dVal[2] );*/
#endif
}
/* fprintf( stderr, "\tCalculated actual bond orders ...\n" );
fprintf( stderr, "j\tDelta\tDelta_e\tDelta_boc\tnlp"
"\tDelta_lp\tClp\tdDelta_lp\n" );
fprintf( stderr, "Atom\tDelta\t\tDelta_e\t\tDelta_boc\tnlp"
"\t\tnlp_opt\t\tDelta_lp\tClp\t\tdDelta_lp\n" );*/
p_lp1 = system->reax_param.gp.l[15];
/* Calculate some helper variables that are used at many places
throughout force calculations */
for( j = 0; j < system->N; ++j ){
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
sbp_j = &(system->reax_param.sbp[ type_j ]);
workspace->Delta[j] = workspace->total_bond_order[j] - sbp_j->valency;
@ -1155,11 +582,7 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
workspace->nlp[j] = explp1 - (int)(workspace->Delta_e[j] / 2.0);
workspace->Delta_lp[j] = sbp_j->nlp_opt - workspace->nlp[j];
workspace->Clp[j] = 2.0 * p_lp1 * explp1 * (2.0 + workspace->vlpex[j]);
/* Adri uses different dDelta_lp values than the ones in notes... */
workspace->dDelta_lp[j] = workspace->Clp[j];
//workspace->dDelta_lp[j] = workspace->Clp[j] + (0.5-workspace->Clp[j]) *
//((fabs(workspace->Delta_e[j]/2.0 -
// (int)(workspace->Delta_e[j]/2.0)) < 0.1) ? 1 : 0 );
if( sbp_j->mass > 21.0 ) {
workspace->nlp_temp[j] = 0.5 * (sbp_j->valency_e - sbp_j->valency);
@ -1171,49 +594,7 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
workspace->Delta_lp_temp[j] = sbp_j->nlp_opt - workspace->nlp_temp[j];
workspace->dDelta_lp_temp[j] = workspace->Clp[j];
}
}
#if defined(TEST_ENERGIES) || defined(TEST_FORCES)
Print_Bond_List( system, control, data, lists, out_control);
#endif
}
#if defined(DEPRECATED)
/* Locate j on i's list.
This function assumes that j is there for sure!
And this is the case given our method of neighbor generation*/
int Locate_Symmetric_Bond( reax_list *bonds, int i, int j )
{
int start = Start_Index(i, bonds);
int end = End_Index(i, bonds);
int mid = (start + end) / 2;
int mid_nbr;
while( (mid_nbr = bonds->select.bond_list[mid].nbr) != j ) {
/*fprintf( stderr, "\tstart: %d end: %d mid: %d\n",
start, end, mid );*/
if( mid_nbr < j )
start = mid+1;
else end = mid - 1;
mid = (start + end) / 2;
}
return mid;
}
inline void Copy_Bond_Order_Data( bond_order_data *dest, bond_order_data *src )
{
dest->BO = src->BO;
dest->BO_s = src->BO_s;
dest->BO_pi = src->BO_pi;
dest->BO_pi2 = src->BO_pi2;
rvec_Scale( dest->dBOp, -1.0, src->dBOp );
rvec_Scale( dest->dln_BOp_s, -1.0, src->dln_BOp_s );
rvec_Scale( dest->dln_BOp_pi, -1.0, src->dln_BOp_pi );
rvec_Scale( dest->dln_BOp_pi2, -1.0, src->dln_BOp_pi2 );
}
#endif

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

@ -25,20 +25,11 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "bonds.h"
#include "bond_orders.h"
#include "list.h"
#include "tool_box.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_bonds.h"
#include "reaxc_bond_orders.h"
#include "reaxc_list.h"
#include "reaxc_tool_box.h"
#include "reaxc_vector.h"
#endif
void Bonds( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
@ -101,19 +92,6 @@ void Bonds( reax_system *system, control_params *control,
bo_ij->Cdbopi -= (CEbo + twbp->De_p);
bo_ij->Cdbopi2 -= (CEbo + twbp->De_pp);
#ifdef TEST_ENERGY
//fprintf( out_control->ebond, "%6d%6d%24.15e%24.15e%24.15e\n",
fprintf( out_control->ebond, "%6d%6d%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
bo_ij->BO, ebond, data->my_en.e_bond );
#endif
#ifdef TEST_FORCES
Add_dBO( system, lists, i, pj, CEbo, workspace->f_be );
Add_dBOpinpi2( system, lists, i, pj,
-(CEbo + twbp->De_p), -(CEbo + twbp->De_pp),
workspace->f_be, workspace->f_be );
#endif
/* Stabilisation terminal triple bond */
if( bo_ij->BO >= 1.00 ) {
if( gp37 == 2 ||
@ -142,17 +120,6 @@ void Bonds( reax_system *system, control_params *control,
bo_ij->Cdbo += decobdbo;
workspace->CdDelta[i] += decobdboua;
workspace->CdDelta[j] += decobdboub;
#ifdef TEST_ENERGY
//fprintf( out_control->ebond,
// "%6d%6d%24.15e%24.15e%24.15e%24.15e\n",
// system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
// estriph, decobdbo, decobdboua, decobdboub );
#endif
#ifdef TEST_FORCES
Add_dBO( system, lists, i, pj, decobdbo, workspace->f_be );
Add_dDelta( system, lists, i, decobdboua, workspace->f_be );
Add_dDelta( system, lists, j, decobdboub, workspace->f_be );
#endif
}
}
}

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

@ -25,14 +25,8 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "control.h"
#include "tool_box.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_control.h"
#include "reaxc_tool_box.h"
#endif
char Read_Control_File( char *control_file, control_params* control,
output_controls *out_control )
@ -122,7 +116,6 @@ char Read_Control_File( char *control_file, control_params* control,
while (!feof(fp)) {
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
//fprintf( stderr, "%s\n", s );
if( strcmp(tmp[0], "simulation_name") == 0 ) {
strcpy( control->sim_name, tmp[1] );
@ -152,13 +145,6 @@ char Read_Control_File( char *control_file, control_params* control,
control->nprocs = control->procs_by_dim[0]*control->procs_by_dim[1]*
control->procs_by_dim[2];
}
//else if( strcmp(tmp[0], "restart") == 0 ) {
// ival = atoi(tmp[1]);
// control->restart = ival;
//}
//else if( strcmp(tmp[0], "restart_from") == 0 ) {
// strcpy( control->restart_from, tmp[1] );
//}
else if( strcmp(tmp[0], "random_vel") == 0 ) {
ival = atoi(tmp[1]);
control->random_vel = ival;
@ -395,13 +381,5 @@ char Read_Control_File( char *control_file, control_params* control,
fclose(fp);
// fprintf( stderr,"%d %d %10.5f %d %10.5f %10.5f\n",
// control->ensemble, control->nsteps, control->dt,
// control->tabulate, control->T, control->P );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "control file read\n" );
#endif
return SUCCESS;
}

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

@ -117,7 +117,7 @@
#define MAX_BOND 20
#define MAXREAXBOND 24 /* used in fix_reaxc_bonds.cpp and pair_reax_c.cpp */
#define MAXSPECBOND 24 /* used in fix_species.cpp and pair_reax_c.cpp */
#define MAXSPECBOND 24 /* used in fix_reaxc_species.cpp and pair_reax_c.cpp */
/******************* ENUMERATIONS *************************/
enum geo_formats { CUSTOM, PDB, ASCII_RESTART, BINARY_RESTART, GF_N };

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

@ -26,14 +26,8 @@
#include "pair_reax_c.h"
#include "error.h"
#if defined(PURE_REAX)
#include "ffield.h"
#include "tool_box.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_ffield.h"
#include "reaxc_tool_box.h"
#endif
char Read_Force_Field( char *ffield_file, reax_interaction *reax,
control_params *control )
@ -158,15 +152,8 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
}
}
// vdWaals type: 1: Shielded Morse, no inner-wall
// 2: inner wall, no shielding
// 3: inner wall+shielding
reax->gp.vdw_type = 0;
/* reading single atom parameters */
/* there are 4 or 5 lines of each single atom parameters in ff files,
depending on using lgvdw or not. These parameters later determine
some of the pair and triplet parameters using combination rules. */
for( i = 0; i < reax->num_atom_types; i++ ) {
/* line one */
@ -261,10 +248,6 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
}
else{
reax->gp.vdw_type = 3;
#if defined(DEBUG)
fprintf( stderr, "vdWaals type for element %s: Shielding+inner-wall",
reax->sbp[i].name );
#endif
}
}
else { // No shielding vdWaals parameters present
@ -278,10 +261,6 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
reax->sbp[i].name );
else{
reax->gp.vdw_type = 2;
#if defined(DEBUG)
fprintf( stderr,"vdWaals type for element%s: No Shielding,inner-wall",
reax->sbp[i].name );
#endif
}
}
}
@ -297,10 +276,6 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
reax->sbp[i].name );
else{
reax->gp.vdw_type = 1;
#if defined(DEBUG)
fprintf( stderr,"vdWaals type for element%s: Shielding,no inner-wall",
reax->sbp[i].name );
#endif
}
}
else{
@ -312,10 +287,6 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
}
}
#if defined(DEBUG)
fprintf( stderr, "vdWaals type: %d\n", reax->gp.vdw_type );
#endif
/* Equate vval3 to valf for first-row elements (25/10/2004) */
for( i = 0; i < reax->num_atom_types; i++ )
if( reax->sbp[i].mass < 21 &&
@ -325,7 +296,6 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
reax->sbp[i].valency_val = reax->sbp[i].valency_boc;
}
/* next line is number of two body combination and some comments */
fgets(s,MAX_LINE,fp);
c=Tokenize(s,&tmp);
@ -385,8 +355,6 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
}
}
/* calculating combination rules and filling up remaining fields. */
for (i=0; i < reax->num_atom_types; i++)
for (j=i; j < reax->num_atom_types; j++) {
reax->tbp[i][j].r_s = 0.5 *
@ -486,10 +454,6 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
}
/* next line is number of two body offdiagonal combinations and comments */
/* these are two body offdiagonal terms that are different from the
combination rules defined above */
fgets(s,MAX_LINE,fp);
c=Tokenize(s,&tmp);
l = atoi(tmp[0]);
@ -546,16 +510,11 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
}
}
/* 3-body parameters -
supports multi-well potentials (upto MAX_3BODY_PARAM in mytypes.h) */
/* clear entries first */
for( i = 0; i < reax->num_atom_types; ++i )
for( j = 0; j < reax->num_atom_types; ++j )
for( k = 0; k < reax->num_atom_types; ++k )
reax->thbp[i][j][k].cnt = 0;
/* next line is number of 3-body params and some comments */
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
l = atoi( tmp[0] );
@ -604,15 +563,6 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
}
}
/* 4-body parameters are entered in compact form. i.e. 0-X-Y-0
correspond to any type of pair of atoms in 1 and 4
position. However, explicit X-Y-Z-W takes precedence over the
default description.
supports multi-well potentials (upto MAX_4BODY_PARAM in mytypes.h)
IMPORTANT: for now, directions on how to read multi-entries from ffield
is not clear */
/* clear all entries first */
for( i = 0; i < reax->num_atom_types; ++i )
for( j = 0; j < reax->num_atom_types; ++j )
@ -639,16 +589,11 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
if (j >= 0 && n >= 0) { // this means the entry is not in compact form
if (j < reax->num_atom_types && k < reax->num_atom_types &&
m < reax->num_atom_types && n < reax->num_atom_types) {
/* these flags ensure that this entry take precedence
over the compact form entries */
tor_flag[j][k][m][n] = 1;
tor_flag[n][m][k][j] = 1;
reax->fbp[j][k][m][n].cnt = 1;
reax->fbp[n][m][k][j].cnt = 1;
/* cnt = reax->fbp[j][k][m][n].cnt;
reax->fbp[j][k][m][n].cnt++;
reax->fbp[n][m][k][j].cnt++; */
val = atof(tmp[4]);
reax->fbp[j][k][m][n].prm[0].V1 = val;
@ -677,9 +622,6 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
for( o = 0; o < reax->num_atom_types; o++ ) {
reax->fbp[p][k][m][o].cnt = 1;
reax->fbp[o][m][k][p].cnt = 1;
/* cnt = reax->fbp[p][k][m][o].cnt;
reax->fbp[p][k][m][o].cnt++;
reax->fbp[o][m][k][p].cnt++; */
if (tor_flag[p][k][m][o] == 0) {
reax->fbp[p][k][m][o].prm[0].V1 = atof(tmp[4]);
@ -731,7 +673,6 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
}
}
/* deallocate helper storage */
for( i = 0; i < MAX_TOKENS; i++ )
free( tmp[i] );
@ -755,9 +696,5 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
fclose(fp);
#if defined(DEBUG_FOCUS)
fprintf( stderr, "force field read\n" );
#endif
return SUCCESS;
}

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

@ -25,23 +25,6 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "forces.h"
#include "bond_orders.h"
#include "bonds.h"
#include "basic_comm.h"
#include "hydrogen_bonds.h"
#include "io_tools.h"
#include "list.h"
#include "lookup.h"
#include "multi_body.h"
#include "nonbonded.h"
#include "qEq.h"
#include "tool_box.h"
#include "torsion_angles.h"
#include "valence_angles.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_forces.h"
#include "reaxc_bond_orders.h"
#include "reaxc_bonds.h"
@ -56,7 +39,6 @@
#include "reaxc_torsion_angles.h"
#include "reaxc_valence_angles.h"
#include "reaxc_vector.h"
#endif
interaction_function Interaction_Functions[NUM_INTRS];
@ -91,23 +73,10 @@ void Compute_Bonded_Forces( reax_system *system, control_params *control,
{
int i;
/* Mark beginning of a new timestep in bonded energy files */
#if defined(TEST_ENERGY)
Debug_Marker_Bonded( out_control, data->step );
#endif
/* Implement all force calls as function pointers */
for( i = 0; i < NUM_INTRS; i++ ) {
#if defined(DEBUG)
fprintf( stderr, "p%d: starting f%d\n", system->my_rank, i );
MPI_Barrier( comm );
#endif
(Interaction_Functions[i])( system, control, data, workspace,
lists, out_control );
#if defined(DEBUG)
fprintf( stderr, "p%d: f%d done\n", system->my_rank, i );
MPI_Barrier( comm );
#endif
}
}
@ -117,10 +86,6 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
reax_list **lists, output_controls *out_control,
MPI_Comm comm )
{
/* Mark beginning of a new timestep in nonbonded energy files */
#if defined(TEST_ENERGY)
Debug_Marker_Nonbonded( out_control, data->step );
#endif
/* van der Waals and Coulomb interactions */
if( control->tabulate == 0 )
@ -129,18 +94,9 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
else
Tabulated_vdW_Coulomb_Energy( system, control, data, workspace,
lists, out_control );
#if defined(DEBUG)
fprintf( stderr, "p%d: nonbonded forces done\n", system->my_rank );
MPI_Barrier( comm );
#endif
}
/* this version of Compute_Total_Force computes forces from
coefficients accumulated by all interaction functions.
Saves enormous time & space! */
void Compute_Total_Force( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, mpi_datatypes *mpi_data )
@ -157,33 +113,6 @@ void Compute_Total_Force( reax_system *system, control_params *control,
Add_dBond_to_Forces_NPT( i, pj, data, workspace, lists );
}
//Print_Total_Force( system, data, workspace );
#if defined(PURE_REAX)
/* now all forces are computed to their partially-final values
based on the neighbors information each processor has had.
final values of force on each atom needs to be computed by adding up
all partially-final pieces */
Coll( system, mpi_data, workspace->f, mpi_data->mpi_rvec,
sizeof(rvec)/sizeof(void), rvec_unpacker );
for( i = 0; i < system->n; ++i )
rvec_Copy( system->my_atoms[i].f, workspace->f[i] );
#if defined(TEST_FORCES)
Coll( system, mpi_data, workspace->f_ele, mpi_data->mpi_rvec, rvec_unpacker);
Coll( system, mpi_data, workspace->f_vdw, mpi_data->mpi_rvec, rvec_unpacker);
Coll( system, mpi_data, workspace->f_be, mpi_data->mpi_rvec, rvec_unpacker );
Coll( system, mpi_data, workspace->f_lp, mpi_data->mpi_rvec, rvec_unpacker );
Coll( system, mpi_data, workspace->f_ov, mpi_data->mpi_rvec, rvec_unpacker );
Coll( system, mpi_data, workspace->f_un, mpi_data->mpi_rvec, rvec_unpacker );
Coll( system, mpi_data, workspace->f_ang, mpi_data->mpi_rvec, rvec_unpacker);
Coll( system, mpi_data, workspace->f_coa, mpi_data->mpi_rvec, rvec_unpacker);
Coll( system, mpi_data, workspace->f_pen, mpi_data->mpi_rvec, rvec_unpacker);
Coll( system, mpi_data, workspace->f_hb, mpi_data->mpi_rvec, rvec_unpacker );
Coll( system, mpi_data, workspace->f_tor, mpi_data->mpi_rvec, rvec_unpacker);
Coll( system, mpi_data, workspace->f_con, mpi_data->mpi_rvec, rvec_unpacker);
#endif
#endif
}
void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists,
@ -201,12 +130,8 @@ void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists,
bonds = *lists + BONDS;
for( i = 0; i < N; ++i ) {
// if( i < n ) - we need to update ghost estimates for delayed nbrings
system->my_atoms[i].num_bonds = MAX(Num_Entries(i,bonds)*2, MIN_BONDS);
//if( End_Index(i, bonds) >= Start_Index(i+1, bonds)-2 )
//workspace->realloc.bonds = 1;
if( i < N-1 )
comp = Start_Index(i+1, bonds);
else comp = bonds->num_intrs;
@ -315,7 +240,6 @@ void Init_Forces( reax_system *system, control_params *control,
workspace->bond_mark[i] = 0;
for( i = system->n; i < system->N; ++i ) {
workspace->bond_mark[i] = 1000; // put ghost atoms to an infinite distance
//workspace->done_after[i] = Start_Index( i, far_nbrs );
}
H = workspace->H;
@ -329,6 +253,7 @@ void Init_Forces( reax_system *system, control_params *control,
for( i = 0; i < system->N; ++i ) {
atom_i = &(system->my_atoms[i]);
type_i = atom_i->type;
if (type_i < 0) continue;
start_i = Start_Index(i, far_nbrs);
end_i = End_Index(i, far_nbrs);
btop_i = End_Index( i, bonds );
@ -364,11 +289,6 @@ void Init_Forces( reax_system *system, control_params *control,
nbr_pj = &( far_nbrs->select.far_nbr_list[pj] );
j = nbr_pj->nbr;
atom_j = &(system->my_atoms[j]);
//fprintf( stderr, "%d%d i=%d x_i: %f %f %f,j=%d x_j: %f %f %f, d=%f\n",
// MIN(atom_i->orig_id, atom_j->orig_id),
// MAX(atom_i->orig_id, atom_j->orig_id),
// i, atom_i->x[0], atom_i->x[1], atom_i->x[2],
// j, atom_j->x[0], atom_j->x[1], atom_j->x[2], nbr_pj->d );
if( renbr ) {
if(nbr_pj->d <= cutoff)
flag = 1;
@ -390,6 +310,7 @@ void Init_Forces( reax_system *system, control_params *control,
if( flag ){
type_j = atom_j->type;
if (type_j < 0) continue;
r_ij = nbr_pj->d;
sbp_j = &(system->reax_param.sbp[type_j]);
twbp = &(system->reax_param.tbp[type_i][type_j]);
@ -398,11 +319,6 @@ void Init_Forces( reax_system *system, control_params *control,
/* H matrix entry */
if( j < system->n || atom_i->orig_id < atom_j->orig_id ) {//tryQEq||1
H->entries[Htop].j = j;
//fprintf( stdout, "%d%d %d %d\n",
// MIN(atom_i->orig_id, atom_j->orig_id),
// MAX(atom_i->orig_id, atom_j->orig_id),
// MIN(atom_i->orig_id, atom_j->orig_id),
// MAX(atom_i->orig_id, atom_j->orig_id) );
if( control->tabulate == 0 )
H->entries[Htop].val = Compute_H(r_ij,twbp->gamma,workspace->Tap);
else H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
@ -412,7 +328,6 @@ void Init_Forces( reax_system *system, control_params *control,
/* hydrogen bond lists */
if( control->hbond_cut > 0 && (ihb==1 || ihb==2) &&
nbr_pj->d <= control->hbond_cut ) {
// fprintf( stderr, "%d %d\n", atom1, atom2 );
jhb = sbp_j->p_hbond;
if( ihb == 1 && jhb == 2 ) {
hbonds->select.hbond_list[ihb_top].nbr = j;
@ -444,11 +359,7 @@ void Init_Forces( reax_system *system, control_params *control,
workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
else if( workspace->bond_mark[i] > workspace->bond_mark[j] + 1 ) {
workspace->bond_mark[i] = workspace->bond_mark[j] + 1;
//if( workspace->bond_mark[i] == 1000 )
// workspace->done_after[i] = pj;
}
//fprintf( stdout, "%d%d - %d(%d) %d(%d)\n",
// i , j, i, workspace->bond_mark[i], j, workspace->bond_mark[j] );
}
}
}
@ -461,77 +372,10 @@ void Init_Forces( reax_system *system, control_params *control,
}
}
//fprintf( stderr, "after the first init loop\n" );
/*for( i = system->n; i < system->N; ++i )
if( workspace->bond_mark[i] > 3 ) {
start_i = Start_Index(i, bonds);
end_i = End_Index(i, bonds);
num_bonds -= (end_i - start_i);
Set_End_Index(i, start_i, bonds );
}*/
/*for( i = system->n; i < system->N; ++i ) {
start_i = Start_Index(i, far_nbrs);
end_i = workspace->done_after[i];
if( workspace->bond_mark[i] >= 2 && start_i < end_i ) {
atom_i = &(system->my_atoms[i]);
type_i = atom_i->type;
btop_i = End_Index( i, bonds );
sbp_i = &(system->reax_param.sbp[type_i]);
for( pj = start_i; pj < end_i; ++pj ) {
nbr_pj = &( far_nbrs->select.far_nbr_list[pj] );
j = nbr_pj->nbr;
if( workspace->bond_mark[j] >= 2 && nbr_pj->d <= control->bond_cut ) {
atom_j = &(system->my_atoms[j]);
type_j = atom_j->type;
sbp_j = &(system->reax_param.sbp[type_j]);
twbp = &(system->reax_param.tbp[type_i][type_j]);
if( BOp( workspace, bonds, control->bo_cut,
i , btop_i, nbr_pj, sbp_i, sbp_j, twbp ) ) {
num_bonds += 2;
++btop_i;
if( workspace->bond_mark[j] > workspace->bond_mark[i] + 1 )
workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
else if( workspace->bond_mark[i] > workspace->bond_mark[j] + 1 )
workspace->bond_mark[i] = workspace->bond_mark[j] + 1;
//fprintf( stdout, "%d%d - %d(%d) %d(%d) new\n",
// i , j, i, workspace->bond_mark[i], j, workspace->bond_mark[j] );
}
}
}
Set_End_Index( i, btop_i, bonds );
}
}*/
workspace->realloc.Htop = Htop;
workspace->realloc.num_bonds = num_bonds;
workspace->realloc.num_hbonds = num_hbonds;
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d @ step%d: Htop = %d num_bonds = %d num_hbonds = %d\n",
system->my_rank, data->step, Htop, num_bonds, num_hbonds );
MPI_Barrier( comm );
#endif
#if defined( DEBUG )
Print_Bonds( system, bonds, "debugbonds.out" );
Print_Bond_List2( system, bonds, "pbonds.out" );
Print_Sparse_Matrix( system, H );
for( i = 0; i < H->n; ++i )
for( j = H->start[i]; j < H->end[i]; ++j )
fprintf( stderr, "%d %d %.15e\n",
MIN(system->my_atoms[i].orig_id,
system->my_atoms[H->entries[j].j].orig_id),
MAX(system->my_atoms[i].orig_id,
system->my_atoms[H->entries[j].j].orig_id),
H->entries[j].val );
#endif
Validate_Lists( system, workspace, lists, data->step,
system->n, system->N, system->numH, comm );
}
@ -562,7 +406,6 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
workspace->bond_mark[i] = 0;
for( i = system->n; i < system->N; ++i ) {
workspace->bond_mark[i] = 1000; // put ghost atoms to an infinite distance
//workspace->done_after[i] = Start_Index( i, far_nbrs );
}
num_bonds = 0;
@ -573,6 +416,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
for( i = 0; i < system->N; ++i ) {
atom_i = &(system->my_atoms[i]);
type_i = atom_i->type;
if (type_i < 0) continue;
start_i = Start_Index(i, far_nbrs);
end_i = End_Index(i, far_nbrs);
btop_i = End_Index( i, bonds );
@ -623,6 +467,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
if( flag ) {
type_j = atom_j->type;
if (type_j < 0) continue;
r_ij = nbr_pj->d;
sbp_j = &(system->reax_param.sbp[type_j]);
twbp = &(system->reax_param.tbp[type_i][type_j]);
@ -651,8 +496,6 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
}
}
/* uncorrected bond orders */
if( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
nbr_pj->d <= control->bond_cut &&
BOp( workspace, bonds, control->bo_cut,
@ -664,11 +507,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
else if( workspace->bond_mark[i] > workspace->bond_mark[j] + 1 ) {
workspace->bond_mark[i] = workspace->bond_mark[j] + 1;
//if( workspace->bond_mark[i] == 1000 )
// workspace->done_after[i] = pj;
}
//fprintf( stdout, "%d%d - %d(%d) %d(%d)\n",
// i , j, i, workspace->bond_mark[i], j, workspace->bond_mark[j] );
}
}
}
@ -678,27 +517,10 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
Set_End_Index( atom_i->Hindex, ihb_top, hbonds );
}
/*for( i = system->n; i < system->N; ++i )
if( workspace->bond_mark[i] > 3 ) {
start_i = Start_Index(i, bonds);
end_i = End_Index(i, bonds);
num_bonds -= (end_i - start_i);
Set_End_Index(i, start_i, bonds );
}*/
workspace->realloc.num_bonds = num_bonds;
workspace->realloc.num_hbonds = num_hbonds;
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d @ step%d: num_bonds = %d num_hbonds = %d\n",
system->my_rank, data->step, num_bonds, num_hbonds );
MPI_Barrier( comm );
#endif
#if defined( DEBUG )
Print_Bonds( system, bonds, "debugbonds.out" );
Print_Bond_List2( system, bonds, "pbonds.out" );
#endif
Validate_Lists( system, workspace, lists, data->step,
system->n, system->N, system->numH, comm );
}
@ -736,6 +558,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
for( i = 0; i < system->N; ++i ) {
atom_i = &(system->my_atoms[i]);
type_i = atom_i->type;
if (type_i < 0) continue;
start_i = Start_Index(i, far_nbrs);
end_i = End_Index(i, far_nbrs);
sbp_i = &(system->reax_param.sbp[type_i]);
@ -759,6 +582,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
if(nbr_pj->d <= cutoff) {
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
r_ij = nbr_pj->d;
sbp_j = &(system->reax_param.sbp[type_j]);
twbp = &(system->reax_param.tbp[type_i][type_j]);
@ -818,16 +642,9 @@ void Estimate_Storages( reax_system *system, control_params *control,
for( i = 0; i < system->N; ++i ) {
*num_3body += SQR(bond_top[i]);
//if( i < system->n )
bond_top[i] = MAX( bond_top[i] * 2, MIN_BONDS );
//else bond_top[i] = MAX_BONDS;
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d @ estimate storages: Htop = %d, num_3body = %d\n",
system->my_rank, *Htop, *num_3body );
MPI_Barrier( comm );
#endif
}
@ -838,23 +655,9 @@ void Compute_Forces( reax_system *system, control_params *control,
{
MPI_Comm comm;
int qeq_flag;
#if defined(LOG_PERFORMANCE)
real t_start = 0;
//MPI_Barrier( mpi_data->world );
if( system->my_rank == MASTER_NODE )
t_start = Get_Time( );
#endif
comm = mpi_data->world;
/********* init forces ************/
#if defined(PURE_REAX)
if( control->qeq_freq && (data->step-data->prev_steps)%control->qeq_freq==0 )
qeq_flag = 1;
else qeq_flag = 0;
#elif defined(LAMMPS_REAX)
qeq_flag = 0;
#endif
if( qeq_flag )
Init_Forces( system, control, data, workspace, lists, out_control, comm );
@ -862,79 +665,16 @@ void Compute_Forces( reax_system *system, control_params *control,
Init_Forces_noQEq( system, control, data, workspace,
lists, out_control, comm );
#if defined(LOG_PERFORMANCE)
//MPI_Barrier( mpi_data->world );
if( system->my_rank == MASTER_NODE )
Update_Timing_Info( &t_start, &(data->timing.init_forces) );
#endif
/********* bonded interactions ************/
Compute_Bonded_Forces( system, control, data, workspace,
lists, out_control, mpi_data->world );
#if defined(LOG_PERFORMANCE)
//MPI_Barrier( mpi_data->world );
if( system->my_rank == MASTER_NODE )
Update_Timing_Info( &t_start, &(data->timing.bonded) );
#endif
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d @ step%d: completed bonded\n",
system->my_rank, data->step );
MPI_Barrier( mpi_data->world );
#endif
/**************** qeq ************************/
#if defined(PURE_REAX)
if( qeq_flag )
QEq( system, control, data, workspace, out_control, mpi_data );
#if defined(LOG_PERFORMANCE)
//MPI_Barrier( mpi_data->world );
if( system->my_rank == MASTER_NODE )
Update_Timing_Info( &t_start, &(data->timing.qEq) );
#endif
#if defined(DEBUG_FOCUS)
fprintf(stderr, "p%d @ step%d: qeq completed\n", system->my_rank, data->step);
MPI_Barrier( mpi_data->world );
#endif
#endif //PURE_REAX
/********* nonbonded interactions ************/
Compute_NonBonded_Forces( system, control, data, workspace,
lists, out_control, mpi_data->world );
#if defined(LOG_PERFORMANCE)
//MPI_Barrier( mpi_data->world );
if( system->my_rank == MASTER_NODE )
Update_Timing_Info( &t_start, &(data->timing.nonb) );
#endif
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d @ step%d: nonbonded forces completed\n",
system->my_rank, data->step );
MPI_Barrier( mpi_data->world );
#endif
/*********** total force ***************/
Compute_Total_Force( system, control, data, workspace, lists, mpi_data );
#if defined(LOG_PERFORMANCE)
//MPI_Barrier( mpi_data->world );
if( system->my_rank == MASTER_NODE )
Update_Timing_Info( &t_start, &(data->timing.bonded) );
#endif
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d @ step%d: total forces computed\n",
system->my_rank, data->step );
//Print_Total_Force( system, data, workspace );
MPI_Barrier( mpi_data->world );
#endif
#if defined(TEST_FORCES)
Print_Force_Files( system, control, data, workspace,
lists, out_control, mpi_data );
#endif
}

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

@ -25,20 +25,11 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "hydrogen_bonds.h"
#include "bond_orders.h"
#include "list.h"
#include "valence_angles.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_hydrogen_bonds.h"
#include "reaxc_bond_orders.h"
#include "reaxc_list.h"
#include "reaxc_valence_angles.h"
#include "reaxc_vector.h"
#endif
void Hydrogen_Bonds( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
@ -55,7 +46,6 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
real 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;
// rtensor temp_rtensor, total_rtensor;
hbond_parameters *hbp;
bond_order_data *bo_ij;
bond_data *pbond_ij;
@ -72,15 +62,9 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
hbonds = (*lists) + HBONDS;
hbond_list = hbonds->select.hbond_list;
/* loops below discover the Hydrogen bonds between i-j-k triplets.
here j is H atom and there has to be some bond between i and j.
Hydrogen bond is between j and k.
so in this function i->X, j->H, k->Z when we map
variables onto the ones in the handout.*/
for( j = 0; j < system->n; ++j )
/* j has to be of type H */
if (system->my_atoms[j].type < 0) continue;
if( system->reax_param.sbp[system->my_atoms[j].type].p_hbond == 1 ) {
/*set j's variables */
type_j = system->my_atoms[j].type;
start_j = Start_Index(j, bonds);
end_j = End_Index(j, bonds);
@ -91,21 +75,20 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
for( pi = start_j; pi < end_j; ++pi ) {
pbond_ij = &( bond_list[pi] );
i = pbond_ij->nbr;
bo_ij = &(pbond_ij->bo_data);
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
bo_ij = &(pbond_ij->bo_data);
if( system->reax_param.sbp[type_i].p_hbond == 2 &&
bo_ij->BO >= HB_THRESHOLD )
hblist[top++] = pi;
}
// fprintf( stderr, "j: %d, top: %d, hb_start_j: %d, hb_end_j:%d\n",
// j, top, hb_start_j, hb_end_j );
for( pk = hb_start_j; pk < hb_end_j; ++pk ) {
/* set k's varibles */
k = hbond_list[pk].nbr;
type_k = system->my_atoms[k].type;
if (type_k < 0) continue;
nbr_jk = hbond_list[pk].ptr;
r_jk = nbr_jk->d;
rvec_Scale( dvec_jk, hbond_list[pk].scl, nbr_jk->dvec );
@ -118,6 +101,7 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
if( system->my_atoms[i].orig_id != system->my_atoms[k].orig_id ) {
bo_ij = &(pbond_ij->bo_data);
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
r_ij = pbond_ij->d;
hbp = &(system->reax_param.hbp[ type_i ][ type_j ][ type_k ]);
++num_hb_intrs;
@ -146,12 +130,6 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
CEhb3 = -hbp->p_hb3 *
(-hbp->r0_hb / SQR(r_jk) + 1.0 / hbp->r0_hb) * e_hb;
/*fprintf( stdout,
"%6d%6d%6d%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n",
system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
r_jk, theta, hbp->p_hb1, exp_hb2, hbp->p_hb3, hbp->r0_hb,
exp_hb3, sin_xhz4, e_hb ); */
/* hydrogen bond forces */
bo_ij->Cdbo += CEhb1; // dbo term
@ -166,8 +144,6 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
rvec_ScaledAdd( workspace->f[k], +CEhb3/r_jk, dvec_jk );
}
else {
/* for pressure coupling, terms that are not related to bond order
derivatives are added directly into pressure vector/tensor */
rvec_Scale( force, +CEhb2, dcos_theta_di ); // dcos terms
rvec_Add( workspace->f[i], force );
rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
@ -202,41 +178,8 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
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",
dcos_theta_di[0], dcos_theta_di[1], dcos_theta_di[2],
dcos_theta_dj[0], dcos_theta_dj[1], dcos_theta_dj[2],
dcos_theta_dk[0], dcos_theta_dk[1], dcos_theta_dk[2]);
fprintf( out_control->ehb, "%24.15e%24.15e%24.15e\n",
CEhb1, CEhb2, CEhb3 ); */
fprintf( out_control->ehb,
//"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
r_jk, theta, bo_ij->BO, e_hb, data->my_en.e_hb );
#endif
#ifdef TEST_FORCES
Add_dBO( system, lists, j, pi, +CEhb1, workspace->f_hb ); //dbo term
// dcos terms
rvec_ScaledAdd( workspace->f_hb[i], +CEhb2, dcos_theta_di );
rvec_ScaledAdd( workspace->f_hb[j], +CEhb2, dcos_theta_dj );
rvec_ScaledAdd( workspace->f_hb[k], +CEhb2, dcos_theta_dk );
// dr terms
rvec_ScaledAdd( workspace->f_hb[j], -CEhb3/r_jk, dvec_jk );
rvec_ScaledAdd( workspace->f_hb[k], +CEhb3/r_jk, dvec_jk );
#endif
}
}
}
}
#if defined(DEBUG)
fprintf( stderr, "Number of hydrogen bonds: %d\n", num_hb_intrs );
fprintf( stderr, "Hydrogen Bond Energy: %g\n", data->my_en.e_hb );
fprintf( stderr, "hydbonds: ext_press (%24.15e %24.15e %24.15e)\n",
data->ext_press[0], data->ext_press[1], data->ext_press[2] );
#endif
}

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

@ -25,24 +25,6 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "init_md.h"
#include "allocate.h"
#include "box.h"
#include "comm_tools.h"
#include "forces.h"
#include "grid.h"
#include "integrate.h"
#include "io_tools.h"
#include "list.h"
#include "lookup.h"
#include "neighbors.h"
#include "random.h"
#include "reset_tools.h"
#include "system_props.h"
#include "tool_box.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_init_md.h"
#include "reaxc_allocate.h"
#include "reaxc_forces.h"
@ -53,228 +35,7 @@
#include "reaxc_system_props.h"
#include "reaxc_tool_box.h"
#include "reaxc_vector.h"
#endif
#if defined(PURE_REAX)
/************************ initialize system ************************/
int Reposition_Atoms( reax_system *system, control_params *control,
simulation_data *data, mpi_datatypes *mpi_data,
char *msg )
{
int i;
rvec dx;
/* reposition atoms */
if( control->reposition_atoms == 0 ) { //fit atoms to periodic box
rvec_MakeZero( dx );
}
else if( control->reposition_atoms == 1 ) { //put center of mass to center
rvec_Scale( dx, 0.5, system->big_box.box_norms );
rvec_ScaledAdd( dx, -1., data->xcm );
}
else if( control->reposition_atoms == 2 ) { //put center of mass to origin
rvec_Scale( dx, -1., data->xcm );
}
else {
strcpy( msg, "reposition_atoms: invalid option" );
return FAILURE;
}
for( i = 0; i < system->n; ++i )
// Inc_on_T3_Gen( system->my_atoms[i].x, dx, &(system->big_box) );
rvec_Add( system->my_atoms[i].x, dx );
return SUCCESS;
}
void Generate_Initial_Velocities( reax_system *system, real T )
{
int i;
real m, scale, norm;
if( T <= 0.1 ) {
for( i = 0; i < system->n; i++ )
rvec_MakeZero( system->my_atoms[i].v );
}
else {
Randomize();
for( i = 0; i < system->n; i++ ) {
rvec_Random( system->my_atoms[i].v );
norm = rvec_Norm_Sqr( system->my_atoms[i].v );
m = system->reax_param.sbp[ system->my_atoms[i].type ].mass;
scale = sqrt( m * norm / (3.0 * K_B * T) );
rvec_Scale( system->my_atoms[i].v, 1./scale, system->my_atoms[i].v );
// fprintf( stderr, "v = %f %f %f\n",
// system->my_atoms[i].v[0],
// system->my_atoms[i].v[1],
// system->my_atoms[i].v[2] );
// fprintf( stderr, "scale = %f\n", scale );
// fprintf( stderr, "v = %f %f %f\n",
// system->my_atoms[i].v[0],
// system->my_atoms[i].v[1],
// system->my_atoms[i].v[2] );
}
}
}
int Init_System( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
mpi_datatypes *mpi_data, char *msg )
{
int i;
reax_atom *atom;
int nrecv[MAX_NBRS];
Setup_New_Grid( system, control, mpi_data->world );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d GRID:\n", system->my_rank );
Print_Grid( &(system->my_grid), stderr );
#endif
Bin_My_Atoms( system, &(workspace->realloc) );
Reorder_My_Atoms( system, workspace );
/* estimate N and total capacity */
for( i = 0; i < MAX_NBRS; ++i ) nrecv[i] = 0;
system->max_recved = 0;
system->N = SendRecv( system, mpi_data, mpi_data->boundary_atom_type, nrecv,
Estimate_Boundary_Atoms, Unpack_Estimate_Message, 1 );
system->total_cap = MAX( (int)(system->N * SAFE_ZONE), MIN_CAP );
Bin_Boundary_Atoms( system );
/* estimate numH and Hcap */
system->numH = 0;
if( control->hbond_cut > 0 )
for( i = 0; i < system->n; ++i ) {
atom = &(system->my_atoms[i]);
if( system->reax_param.sbp[ atom->type ].p_hbond == 1 )
atom->Hindex = system->numH++;
else atom->Hindex = -1;
}
system->Hcap = MAX( system->numH * SAFER_ZONE, MIN_CAP );
//Allocate_System( system, system->local_cap, system->total_cap, msg );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: n=%d local_cap=%d\n",
system->my_rank, system->n, system->local_cap );
fprintf( stderr, "p%d: N=%d total_cap=%d\n",
system->my_rank, system->N, system->total_cap );
fprintf( stderr, "p%d: numH=%d H_cap=%d\n",
system->my_rank, system->numH, system->Hcap );
#endif
// if( Reposition_Atoms( system, control, data, mpi_data, msg ) == FAILURE )
// return FAILURE;
/* initialize velocities so that desired init T can be attained */
if( !control->restart || (control->restart && control->random_vel) )
Generate_Initial_Velocities( system, control->T_init );
return SUCCESS;
}
/************************ initialize simulation data ************************/
int Init_Simulation_Data( reax_system *system, control_params *control,
simulation_data *data, mpi_datatypes *mpi_data,
char *msg )
{
Reset_Simulation_Data( data, control->virial );
if( !control->restart )
data->step = data->prev_steps = 0;
Compute_Total_Mass( system, data, mpi_data->comm_mesh3D );
Compute_Center_of_Mass( system, data, mpi_data, mpi_data->comm_mesh3D );
Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
switch( control->ensemble ){
case NVE:
data->N_f = 3 * system->bigN;
Evolve = Velocity_Verlet_NVE;
break;
case bNVT:
data->N_f = 3 * system->bigN + 1;
Evolve = Velocity_Verlet_Berendsen_NVT;
break;
case nhNVT:
fprintf( stderr, "WARNING: Nose-Hoover NVT is still under testing.\n" );
//return FAILURE;
data->N_f = 3 * system->bigN + 1;
Evolve = Velocity_Verlet_Nose_Hoover_NVT_Klein;
if( !control->restart || (control->restart && control->random_vel) ) {
data->therm.G_xi = control->Tau_T *
(2.0 * data->sys_en.e_kin - data->N_f * K_B * control->T );
data->therm.v_xi = data->therm.G_xi * control->dt;
data->therm.v_xi_old = 0;
data->therm.xi = 0;
}
break;
case sNPT: /* Semi-Isotropic NPT */
data->N_f = 3 * system->bigN + 4;
Evolve = Velocity_Verlet_Berendsen_NPT;
if( !control->restart )
Reset_Pressures( data );
break;
case iNPT: /* Isotropic NPT */
data->N_f = 3 * system->bigN + 2;
Evolve = Velocity_Verlet_Berendsen_NPT;
if( !control->restart )
Reset_Pressures( data );
break;
case NPT: /* Anisotropic NPT */
strcpy( msg, "init_simulation_data: option not yet implemented" );
return FAILURE;
data->N_f = 3 * system->bigN + 9;
Evolve = Velocity_Verlet_Berendsen_NPT;
/*if( !control->restart ) {
data->therm.G_xi = control->Tau_T *
(2.0 * data->my_en.e_Kin - data->N_f * K_B * control->T );
data->therm.v_xi = data->therm.G_xi * control->dt;
data->iso_bar.eps = 0.33333 * log(system->box.volume);
data->inv_W = 1.0 /
( data->N_f * K_B * control->T * SQR(control->Tau_P) );
Compute_Pressure( system, control, data, out_control );
}*/
break;
default:
strcpy( msg, "init_simulation_data: ensemble not recognized" );
return FAILURE;
}
/* initialize the timer(s) */
MPI_Barrier( mpi_data->world ); // wait for everyone to come here
if( system->my_rank == MASTER_NODE ) {
data->timing.start = Get_Time( );
#if defined(LOG_PERFORMANCE)
Reset_Timing( &data->timing );
#endif
}
#if defined(DEBUG)
fprintf( stderr, "data->N_f: %8.3f\n", data->N_f );
#endif
return SUCCESS;
}
#elif defined(LAMMPS_REAX)
int Init_System( reax_system *system, control_params *control, char *msg )
{
int i;
@ -294,21 +55,13 @@ int Init_System( reax_system *system, control_params *control, char *msg )
if( control->hbond_cut > 0 )
for( i = 0; i < system->n; ++i ) {
atom = &(system->my_atoms[i]);
if (atom->type < 0) continue;
if( system->reax_param.sbp[ atom->type ].p_hbond == 1 )
atom->Hindex = system->numH++;
else atom->Hindex = -1;
}
system->Hcap = (int)(MAX( system->numH * saferzone, mincap ));
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: n=%d local_cap=%d\n",
system->my_rank, system->n, system->local_cap );
fprintf( stderr, "p%d: N=%d total_cap=%d\n",
system->my_rank, system->N, system->total_cap );
fprintf( stderr, "p%d: numH=%d H_cap=%d\n",
system->my_rank, system->numH, system->Hcap );
#endif
return SUCCESS;
}
@ -321,22 +74,13 @@ int Init_Simulation_Data( reax_system *system, control_params *control,
/* initialize the timer(s) */
if( system->my_rank == MASTER_NODE ) {
data->timing.start = Get_Time( );
#if defined(LOG_PERFORMANCE)
Reset_Timing( &data->timing );
#endif
}
//if( !control->restart )
data->step = data->prev_steps = 0;
return SUCCESS;
}
#endif
/************************ initialize workspace ************************/
/* Initialize Taper params */
void Init_Taper( control_params *control, storage *workspace, MPI_Comm comm )
{
real d1, d7;
@ -399,241 +143,14 @@ int Init_Workspace( reax_system *system, control_params *control,
int Init_MPI_Datatypes( reax_system *system, storage *workspace,
mpi_datatypes *mpi_data, MPI_Comm comm, char *msg )
{
#if defined(PURE_REAX)
int i, block[11];
MPI_Aint base, disp[11];
MPI_Datatype type[11];
mpi_atom sample;
boundary_atom b_sample;
restart_atom r_sample;
rvec rvec_sample;
rvec2 rvec2_sample;
#endif
/* setup the world */
mpi_data->world = comm;
MPI_Comm_size( comm, &(system->wsize) );
#if defined(PURE_REAX)
/* init mpi buffers */
mpi_data->in1_buffer = NULL;
mpi_data->in2_buffer = NULL;
/* mpi_atom - [orig_id, imprt_id, type, num_bonds, num_hbonds, name,
x, v, f_old, s, t] */
block[0] = block[1] = block[2] = block[3] = block[4] = 1;
block[5] = 8;
block[6] = block[7] = block[8] = 3;
block[9] = block[10] = 4;
MPI_Address( &(sample.orig_id), disp + 0 );
MPI_Address( &(sample.imprt_id), disp + 1 );
MPI_Address( &(sample.type), disp + 2 );
MPI_Address( &(sample.num_bonds), disp + 3 );
MPI_Address( &(sample.num_hbonds), disp + 4 );
MPI_Address( &(sample.name), disp + 5 );
MPI_Address( &(sample.x[0]), disp + 6 );
MPI_Address( &(sample.v[0]), disp + 7 );
MPI_Address( &(sample.f_old[0]), disp + 8 );
MPI_Address( &(sample.s[0]), disp + 9 );
MPI_Address( &(sample.t[0]), disp + 10 );
base = (MPI_Aint)(&(sample));
for( i = 0; i < 11; ++i ) disp[i] -= base;
type[0] = type[1] = type[2] = type[3] = type[4] = MPI_INT;
type[5] = MPI_CHAR;
type[6] = type[7] = type[8] = type[9] = type[10] = MPI_DOUBLE;
MPI_Type_struct( 11, block, disp, type, &(mpi_data->mpi_atom_type) );
MPI_Type_commit( &(mpi_data->mpi_atom_type) );
/* boundary_atom - [orig_id, imprt_id, type, num_bonds, num_hbonds, x] */
block[0] = block[1] = block[2] = block[3] = block[4] = 1;
block[5] = 3;
MPI_Address( &(b_sample.orig_id), disp + 0 );
MPI_Address( &(b_sample.imprt_id), disp + 1 );
MPI_Address( &(b_sample.type), disp + 2 );
MPI_Address( &(b_sample.num_bonds), disp + 3 );
MPI_Address( &(b_sample.num_hbonds), disp + 4 );
MPI_Address( &(b_sample.x[0]), disp + 5 );
base = (MPI_Aint)(&(b_sample));
for( i = 0; i < 6; ++i ) disp[i] -= base;
type[0] = type[1] = type[2] = type[3] = type[4] = MPI_INT;
type[5] = MPI_DOUBLE;
MPI_Type_struct( 6, block, disp, type, &(mpi_data->boundary_atom_type) );
MPI_Type_commit( &(mpi_data->boundary_atom_type) );
/* mpi_rvec */
block[0] = 3;
MPI_Address( &(rvec_sample[0]), disp + 0 );
base = disp[0];
for( i = 0; i < 1; ++i ) disp[i] -= base;
type[0] = MPI_DOUBLE;
MPI_Type_struct( 1, block, disp, type, &(mpi_data->mpi_rvec) );
MPI_Type_commit( &(mpi_data->mpi_rvec) );
/* mpi_rvec2 */
block[0] = 2;
MPI_Address( &(rvec2_sample[0]), disp + 0 );
base = disp[0];
for( i = 0; i < 1; ++i ) disp[i] -= base;
type[0] = MPI_DOUBLE;
MPI_Type_struct( 1, block, disp, type, &(mpi_data->mpi_rvec2) );
MPI_Type_commit( &(mpi_data->mpi_rvec2) );
/* restart_atom - [orig_id, type, name[8], x, v] */
block[0] = block[1] = 1 ;
block[2] = 8;
block[3] = block[4] = 3;
MPI_Address( &(r_sample.orig_id), disp + 0 );
MPI_Address( &(r_sample.type), disp + 1 );
MPI_Address( &(r_sample.name), disp + 2 );
MPI_Address( &(r_sample.x[0]), disp + 3 );
MPI_Address( &(r_sample.v[0]), disp + 4 );
base = (MPI_Aint)(&(r_sample));
for( i = 0; i < 5; ++i ) disp[i] -= base;
type[0] = type[1] = MPI_INT;
type[2] = MPI_CHAR;
type[3] = type[4] = MPI_DOUBLE;
MPI_Type_struct( 5, block, disp, type, &(mpi_data->restart_atom_type) );
MPI_Type_commit( &(mpi_data->restart_atom_type) );
#endif
return SUCCESS;
}
/********************** allocate lists *************************/
#if defined(PURE_REAX)
int Init_Lists( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
mpi_datatypes *mpi_data, char *msg )
{
int i, num_nbrs;
int total_hbonds, total_bonds, bond_cap, num_3body, cap_3body, Htop;
int *hb_top, *bond_top;
MPI_Comm comm;
comm = mpi_data->world;
//for( i = 0; i < MAX_NBRS; ++i ) nrecv[i] = system->my_nbrs[i].est_recv;
//system->N = SendRecv( system, mpi_data, mpi_data->boundary_atom_type, nrecv,
// Sort_Boundary_Atoms, Unpack_Exchange_Message, 1 );
num_nbrs = Estimate_NumNeighbors( system, lists );
if(!Make_List( system->total_cap, num_nbrs, TYP_FAR_NEIGHBOR,
*lists+FAR_NBRS, comm )){
fprintf(stderr, "Problem in initializing far nbrs list. Terminating!\n");
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: allocated far_nbrs: num_far=%d, space=%dMB\n",
system->my_rank, num_nbrs,
(int)(num_nbrs*sizeof(far_neighbor_data)/(1024*1024)) );
#endif
Generate_Neighbor_Lists( system, data, workspace, lists );
bond_top = (int*) calloc( system->total_cap, sizeof(int) );
hb_top = (int*) calloc( system->local_cap, sizeof(int) );
Estimate_Storages( system, control, lists,
&Htop, hb_top, bond_top, &num_3body, comm );
Allocate_Matrix( &(workspace->H), system->local_cap, Htop, comm );
workspace->L = NULL;
workspace->U = NULL;
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: allocated H matrix: Htop=%d, space=%dMB\n",
system->my_rank, Htop,
(int)(Htop * sizeof(sparse_matrix_entry) / (1024*1024)) );
#endif
if( control->hbond_cut > 0 ) {
/* init H indexes */
total_hbonds = 0;
for( i = 0; i < system->n; ++i ) {
system->my_atoms[i].num_hbonds = hb_top[i];
total_hbonds += hb_top[i];
}
total_hbonds = MAX( total_hbonds*SAFER_ZONE, MIN_CAP*MIN_HBONDS );
if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
*lists+HBONDS, comm ) ) {
fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: allocated hbonds: total_hbonds=%d, space=%dMB\n",
system->my_rank, total_hbonds,
(int)(total_hbonds*sizeof(hbond_data)/(1024*1024)) );
#endif
}
/* bonds list */
//Allocate_Bond_List( system->N, bond_top, (*lists)+BONDS );
//num_bonds = bond_top[system->N-1];
total_bonds = 0;
for( i = 0; i < system->N; ++i ) {
system->my_atoms[i].num_bonds = bond_top[i];
total_bonds += bond_top[i];
}
bond_cap = MAX( total_bonds*SAFE_ZONE, MIN_CAP*MIN_BONDS );
if( !Make_List( system->total_cap, bond_cap, TYP_BOND,
*lists+BONDS, comm ) ) {
fprintf( stderr, "not enough space for bonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: allocated bonds: total_bonds=%d, space=%dMB\n",
system->my_rank, bond_cap,
(int)(bond_cap*sizeof(bond_data)/(1024*1024)) );
#endif
/* 3bodies list */
cap_3body = MAX( num_3body*SAFE_ZONE, MIN_3BODIES );
if( !Make_List( bond_cap, cap_3body, TYP_THREE_BODY,
*lists+THREE_BODIES, comm ) ){
fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: allocated 3-body list: num_3body=%d, space=%dMB\n",
system->my_rank, cap_3body,
(int)(cap_3body*sizeof(three_body_interaction_data)/(1024*1024)) );
#endif
#if defined(TEST_FORCES)
if( !Make_List( system->total_cap, bond_cap*8, TYP_DDELTA,
*lists+DDELTAS, comm ) ) {
fprintf( stderr, "Problem in initializing dDelta list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
fprintf( stderr, "p%d: allocated dDelta list: num_ddelta=%d space=%ldMB\n",
system->my_rank, bond_cap*30,
bond_cap*8*sizeof(dDelta_data)/(1024*1024) );
if( !Make_List( bond_cap, bond_cap*50, TYP_DBO, *lists+DBOS, comm ) ) {
fprintf( stderr, "Problem in initializing dBO list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
fprintf( stderr, "p%d: allocated dbond list: num_dbonds=%d space=%ldMB\n",
system->my_rank, bond_cap*MAX_BONDS*3,
bond_cap*MAX_BONDS*3*sizeof(dbond_data)/(1024*1024) );
#endif
free( hb_top );
free( bond_top );
return SUCCESS;
}
#elif defined(LAMMPS_REAX)
int Init_Lists( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
mpi_datatypes *mpi_data, char *msg )
@ -668,16 +185,8 @@ int Init_Lists( reax_system *system, control_params *control,
fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: allocated hbonds: total_hbonds=%d, space=%dMB\n",
system->my_rank, total_hbonds,
(int)(total_hbonds*sizeof(hbond_data)/(1024*1024)) );
#endif
}
/* bonds list */
//Allocate_Bond_List( system->N, bond_top, (*lists)+BONDS );
//num_bonds = bond_top[system->N-1];
total_bonds = 0;
for( i = 0; i < system->N; ++i ) {
system->my_atoms[i].num_bonds = bond_top[i];
@ -690,11 +199,6 @@ int Init_Lists( reax_system *system, control_params *control,
fprintf( stderr, "not enough space for bonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: allocated bonds: total_bonds=%d, space=%dMB\n",
system->my_rank, bond_cap,
(int)(bond_cap*sizeof(bond_data)/(1024*1024)) );
#endif
/* 3bodies list */
cap_3body = (int)(MAX( num_3body*safezone, MIN_3BODIES ));
@ -703,136 +207,13 @@ int Init_Lists( reax_system *system, control_params *control,
fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: allocated 3-body list: num_3body=%d, space=%dMB\n",
system->my_rank, cap_3body,
(int)(cap_3body*sizeof(three_body_interaction_data)/(1024*1024)) );
#endif
#if defined(TEST_FORCES)
if( !Make_List( system->total_cap, bond_cap*8, TYP_DDELTA,
*lists+DDELTAS, comm ) ) {
fprintf( stderr, "Problem in initializing dDelta list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
fprintf( stderr, "p%d: allocated dDelta list: num_ddelta=%d space=%ldMB\n",
system->my_rank, bond_cap*30,
bond_cap*8*sizeof(dDelta_data)/(1024*1024) );
if( !Make_List( bond_cap, bond_cap*50, TYP_DBO, (*lists)+DBOS, comm ) ) {
fprintf( stderr, "Problem in initializing dBO list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
fprintf( stderr, "p%d: allocated dbond list: num_dbonds=%d space=%ldMB\n",
system->my_rank, bond_cap*MAX_BONDS*3,
bond_cap*MAX_BONDS*3*sizeof(dbond_data)/(1024*1024) );
#endif
free( hb_top );
free( bond_top );
return SUCCESS;
}
#endif
#if defined(PURE_REAX)
void Initialize( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control,
mpi_datatypes *mpi_data )
{
char msg[MAX_STR];
if( Init_MPI_Datatypes( system, workspace, mpi_data, MPI_COMM_WORLD, msg ) ==
FAILURE ) {
fprintf( stderr, "p%d: init_mpi_datatypes: could not create datatypes\n",
system->my_rank );
fprintf( stderr, "p%d: mpi_data couldn't be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized mpi datatypes\n", system->my_rank );
#endif
if( Init_System(system, control, data, workspace, mpi_data, msg) == FAILURE ){
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: system could not be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: system initialized\n", system->my_rank );
#endif
if( Init_Simulation_Data(system, control, data, mpi_data, msg) == FAILURE ) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: sim_data couldn't be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized simulation data\n", system->my_rank );
#endif
if( Init_Workspace( system, control, workspace, mpi_data->world, msg ) ==
FAILURE ) {
fprintf( stderr, "p%d:init_workspace: not enough memory\n",
system->my_rank );
fprintf( stderr, "p%d:workspace couldn't be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized workspace\n", system->my_rank );
#endif
if( Init_Lists( system, control, data, workspace, lists, mpi_data, msg ) ==
FAILURE ) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: system could not be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized lists\n", system->my_rank );
#endif
if(Init_Output_Files(system,control,out_control,mpi_data,msg) == FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: could not open output files! terminating...\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: output files opened\n", system->my_rank );
#endif
if( control->tabulate ) {
if( Init_Lookup_Tables(system,control,workspace,mpi_data,msg) == FAILURE ) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: couldn't create lookup table! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized lookup tables\n", system->my_rank );
#endif
}
Init_Force_Functions( control );
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized force functions\n", system->my_rank );
#endif
/*#ifdef TEST_FORCES
Init_Force_Test_Functions();
fprintf(stderr,"p%d: initialized force test functions\n",system->my_rank);
#endif */
}
#elif defined(LAMMPS_REAX)
void Initialize( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control,
@ -848,9 +229,6 @@ void Initialize( reax_system *system, control_params *control,
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized mpi datatypes\n", system->my_rank );
#endif
if( Init_System(system, control, msg) == FAILURE ){
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
@ -858,9 +236,6 @@ void Initialize( reax_system *system, control_params *control,
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: system initialized\n", system->my_rank );
#endif
if( Init_Simulation_Data( system, control, data, msg ) == FAILURE ) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
@ -868,9 +243,6 @@ void Initialize( reax_system *system, control_params *control,
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized simulation data\n", system->my_rank );
#endif
if( Init_Workspace( system, control, workspace, mpi_data->world, msg ) ==
FAILURE ) {
@ -880,9 +252,6 @@ void Initialize( reax_system *system, control_params *control,
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized workspace\n", system->my_rank );
#endif
if( Init_Lists( system, control, data, workspace, lists, mpi_data, msg ) ==
FAILURE ) {
@ -891,9 +260,6 @@ void Initialize( reax_system *system, control_params *control,
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized lists\n", system->my_rank );
#endif
if( Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
@ -901,9 +267,6 @@ void Initialize( reax_system *system, control_params *control,
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: output files opened\n", system->my_rank );
#endif
if( control->tabulate ) {
if( Init_Lookup_Tables( system, control, workspace, mpi_data, msg ) == FAILURE ) {
@ -912,19 +275,8 @@ void Initialize( reax_system *system, control_params *control,
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
}
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized lookup tables\n", system->my_rank );
#endif
}
Init_Force_Functions( control );
#if defined(DEBUG)
fprintf( stderr, "p%d: initialized force functions\n", system->my_rank );
#endif
/*#if defined(TEST_FORCES)
Init_Force_Test_Functions();
fprintf(stderr,"p%d: initialized force test functions\n",system->my_rank);
#endif*/
}
#endif

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

@ -29,12 +29,6 @@
#include "reaxc_types.h"
#if defined(PURE_REAX)
void Initialize( reax_system*, control_params*, simulation_data*,
storage*, reax_list**, output_controls*, mpi_datatypes* );
#elif defined(LAMMPS_REAX)
void Initialize( reax_system*, control_params*, simulation_data*, storage*,
reax_list**, output_controls*, mpi_datatypes*, MPI_Comm );
#endif
#endif

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

File diff suppressed because it is too large Load Diff

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

@ -25,14 +25,8 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "list.h"
#include "tool_box.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_list.h"
#include "reaxc_tool_box.h"
#endif
/************* allocate list space ******************/
int Make_List(int n, int num_intrs, int type, reax_list *l, MPI_Comm comm)
@ -46,9 +40,6 @@ int Make_List(int n, int num_intrs, int type, reax_list *l, MPI_Comm comm)
l->end_index = (int*) smalloc( n * sizeof(int), "list:end_index", comm );
l->type = type;
#if defined(DEBUG_FOCUS)
fprintf( stderr, "list: n=%d num_intrs=%d type=%d\n", n, num_intrs, type );
#endif
switch(l->type) {
case TYP_VOID:
@ -133,30 +124,3 @@ void Delete_List( reax_list *l, MPI_Comm comm )
}
}
#if defined(PURE_REAX)
inline int Num_Entries( int i, reax_list *l )
{
return l->end_index[i] - l->index[i];
}
inline int Start_Index( int i, reax_list *l )
{
return l->index[i];
}
inline int End_Index( int i, reax_list *l )
{
return l->end_index[i];
}
inline void Set_Start_Index( int i, int val, reax_list *l )
{
l->index[i] = val;
}
inline void Set_End_Index( int i, int val, reax_list *l )
{
l->end_index[i] = val;
}
#endif

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

@ -25,25 +25,17 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "lookup.h"
#include "nonbonded.h"
#include "tool_box.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_lookup.h"
#include "reaxc_nonbonded.h"
#include "reaxc_tool_box.h"
#endif
LR_lookup_table **LR;
/* Fills solution into x. Warning: will modify c and d! */
void Tridiagonal_Solve( const real *a, const real *b,
real *c, real *d, real *x, unsigned int n){
int i;
real id;
/* Modify the coefficients. */
c[0] /= b[0]; /* Division by zero risk. */
d[0] /= b[0]; /* Division by zero would imply a singular matrix. */
for(i = 1; i < n; i++){
@ -52,7 +44,6 @@ void Tridiagonal_Solve( const real *a, const real *b,
d[i] = (d[i] - d[i-1] * a[i])/id;
}
/* Now back substitute. */
x[n - 1] = d[n - 1];
for(i = n - 2; i >= 0; i--)
x[i] = d[i] - c[i] * x[i + 1];
@ -216,17 +207,12 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
fCEclmb = (real*)
smalloc( (control->tabulate+2) * sizeof(real), "lookup:fCEclmb", comm );
/* allocate Long-Range LookUp Table space based on
number of atom types in the ffield file */
LR = (LR_lookup_table**)
scalloc( num_atom_types, sizeof(LR_lookup_table*), "lookup:LR", comm );
for( i = 0; i < num_atom_types; ++i )
LR[i] = (LR_lookup_table*)
scalloc( num_atom_types, sizeof(LR_lookup_table), "lookup:LR[i]", comm );
/* most atom types in ffield file will not exist in the current
simulation. to avoid unnecessary lookup table space, determine
the atom types that exist in the current simulation */
for( i = 0; i < MAX_ATOM_TYPES; ++i )
existing_types[i] = 0;
for( i = 0; i < system->n; ++i )
@ -235,11 +221,8 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
MPI_Allreduce( existing_types, aggregated, MAX_ATOM_TYPES,
MPI_INT, MPI_SUM, mpi_data->world );
/* fill in the lookup table entries for existing atom types.
only lower half should be enough. */
for( i = 0; i < num_atom_types; ++i ) {
if( aggregated[i] ) {
//for( j = 0; j < num_atom_types; ++j )
for( j = i; j < num_atom_types; ++j ) {
if( aggregated[j] ) {
LR[i][j].xmin = 0;

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

@ -25,18 +25,10 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "multi_body.h"
#include "bond_orders.h"
#include "list.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_multi_body.h"
#include "reaxc_bond_orders.h"
#include "reaxc_list.h"
#include "reaxc_vector.h"
#endif
void Atom_Energy( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
@ -51,17 +43,19 @@ void Atom_Energy( reax_system *system, control_params *control,
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_lp1, p_lp2, p_lp3;
real p_ovun2, p_ovun3, p_ovun4, p_ovun5, p_ovun6, p_ovun7, p_ovun8;
real eng_tmp;
real eng_tmp, f_tmp;
int numbonds;
single_body_parameters *sbp_i;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
bond_data *pbond;
bond_order_data *bo_ij;
reax_list *bonds = (*lists) + BONDS;
/* Initialize parameters */
p_lp1 = system->reax_param.gp.l[15];
p_lp3 = system->reax_param.gp.l[5];
p_ovun3 = system->reax_param.gp.l[32];
p_ovun4 = system->reax_param.gp.l[31];
@ -72,6 +66,7 @@ void Atom_Energy( reax_system *system, control_params *control,
for( i = 0; i < system->n; ++i ) {
/* set the parameter pointer */
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
sbp_i = &(system->reax_param.sbp[ type_i ]);
/* lone-pair Energy */
@ -79,37 +74,32 @@ void Atom_Energy( reax_system *system, control_params *control,
expvd2 = exp( -75 * workspace->Delta_lp[i] );
inv_expvd2 = 1. / (1. + expvd2 );
numbonds = 0;
e_lp = 0.0;
for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
numbonds ++;
/* calculate the energy */
data->my_en.e_lp += e_lp =
p_lp2 * workspace->Delta_lp[i] * inv_expvd2;
if (numbonds > 0)
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
if (numbonds > 0) workspace->CdDelta[i] += CElp; // lp - 1st term
/* tally into per-atom energy */
if( system->pair_ptr->evflag)
system->pair_ptr->ev_tally(i,i,system->n,1,e_lp,0.0,0.0,0.0,0.0,0.0);
#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 );
// fprintf( out_control->elp, "%6d%24.15e%24.15e%24.15e\n",
fprintf( out_control->elp, "%6d%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id, workspace->nlp[i],
e_lp, data->my_en.e_lp );
#endif
#ifdef TEST_FORCES
Add_dDelta( system, lists, i, CElp, workspace->f_lp ); // lp - 1st term
#endif
/* correction for C2 */
if( p_lp3 > 0.001 && !strcmp(system->reax_param.sbp[type_i].name, "C") )
for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj ) {
j = bonds->select.bond_list[pj].nbr;
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
if( !strcmp( system->reax_param.sbp[type_j].name, "C" ) ) {
twbp = &( system->reax_param.tbp[type_i][type_j]);
@ -130,15 +120,6 @@ void Atom_Energy( reax_system *system, control_params *control,
if( system->pair_ptr->evflag)
system->pair_ptr->ev_tally(i,j,system->n,1,e_lph,0.0,0.0,0.0,0.0,0.0);
#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,
e_lph, deahu2dbo, deahu2dsbo );
#endif
#ifdef TEST_FORCES
Add_dBO(system, lists, i, pj, deahu2dbo, workspace->f_lp);
Add_dDelta(system, lists, i, deahu2dsbo, workspace->f_lp);
#endif
}
}
}
@ -147,6 +128,7 @@ void Atom_Energy( reax_system *system, control_params *control,
for( i = 0; i < system->n; ++i ) {
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
sbp_i = &(system->reax_param.sbp[ type_i ]);
/* over-coordination energy */
@ -159,18 +141,15 @@ void Atom_Energy( reax_system *system, control_params *control,
for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj ) {
j = bonds->select.bond_list[pj].nbr;
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
bo_ij = &(bonds->select.bond_list[pj].bo_data);
sbp_j = &(system->reax_param.sbp[ type_j ]);
twbp = &(system->reax_param.tbp[ type_i ][ type_j ]);
sum_ovun1 += twbp->p_ovun1 * twbp->De_s * bo_ij->BO;
sum_ovun2 += (workspace->Delta[j] - dfvl*workspace->Delta_lp_temp[j])*
( bo_ij->BO_pi + bo_ij->BO_pi2 );
/*fprintf( stdout, "%4d%4d%12.6f%12.6f%12.6f\n",
i+1, j+1,
dfvl * workspace->Delta_lp_temp[j],
sbp_j->nlp_opt,
workspace->nlp_temp[j] );*/
}
exp_ovun1 = p_ovun3 * exp( p_ovun4 * sum_ovun2 );
@ -205,8 +184,14 @@ void Atom_Energy( reax_system *system, control_params *control,
inv_exp_ovun2n = 1.0 / (1.0 + exp_ovun2n);
inv_exp_ovun8 = 1.0 / (1.0 + exp_ovun8);
data->my_en.e_un += e_un =
-p_ovun5 * (1.0 - exp_ovun6) * inv_exp_ovun2n * inv_exp_ovun8;
numbonds = 0;
e_un = 0.0;
for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
numbonds ++;
if (numbonds > 0)
data->my_en.e_un += e_un =
-p_ovun5 * (1.0 - exp_ovun6) * inv_exp_ovun2n * inv_exp_ovun8;
CEunder1 = inv_exp_ovun2n *
( p_ovun5 * p_ovun6 * exp_ovun6 * inv_exp_ovun8 +
@ -218,18 +203,15 @@ void Atom_Energy( reax_system *system, control_params *control,
/* tally into per-atom energy */
if( system->pair_ptr->evflag) {
eng_tmp = e_ov + e_un;
eng_tmp = e_ov;
if (numbonds > 0) eng_tmp += 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
workspace->CdDelta[i] += CEunder3; // UnCoor - 1st term
#ifdef TEST_FORCES
Add_dDelta( system, lists, i, CEover3, workspace->f_ov ); // OvCoor 2nd
Add_dDelta( system, lists, i, CEunder3, workspace->f_un ); // UnCoor 1st
#endif
if (numbonds > 0) workspace->CdDelta[i] += CEunder3; // UnCoor - 1st term
for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj ) {
pbond = &(bonds->select.bond_list[pj]);
@ -255,72 +237,7 @@ void Atom_Energy( reax_system *system, control_params *control,
bo_ij->Cdbopi2 += CEunder4 *
(workspace->Delta[j] - dfvl*workspace->Delta_lp_temp[j]); // UnCoor-2b
#ifdef TEST_ENERGY
/* fprintf( out_control->eov, "%6d%12.6f\n",
workspace->reverse_map[j],
// CEover1 * twbp->p_ovun1 * twbp->De_s, CEover3,
CEover4 * (1.0 - workspace->dDelta_lp[j]) *
(bo_ij->BO_pi + bo_ij->BO_pi2)
*/// /*CEover4 * (workspace->Delta[j]-workspace->Delta_lp[j])*/);
// fprintf( out_control->eov, "%6d%12.6f\n",
// fprintf( out_control->eov, "%6d%24.15e\n",
// system->my_atoms[j].orig_id,
// CEover1 * twbp->p_ovun1 * twbp->De_s, CEover3,
// CEover4 * (1.0 - workspace->dDelta_lp[j]) *
// (bo_ij->BO_pi + bo_ij->BO_pi2)
// /*CEover4 * (workspace->Delta[j]-workspace->Delta_lp[j])*/);
// CEunder4 * (1.0 - workspace->dDelta_lp[j]) *
// (bo_ij->BO_pi + bo_ij->BO_pi2),
// CEunder4 * (workspace->Delta[j] - workspace->Delta_lp[j]) );
#endif
#ifdef TEST_FORCES
Add_dBO( system, lists, i, pj, CEover1 * twbp->p_ovun1 * twbp->De_s,
workspace->f_ov ); // OvCoor - 1st term
Add_dDelta( system, lists, j,
CEover4 * (1.0 - dfvl*workspace->dDelta_lp[j]) *
(bo_ij->BO_pi + bo_ij->BO_pi2),
workspace->f_ov ); // OvCoor - 3a
Add_dBOpinpi2( system, lists, i, pj,
CEover4 * (workspace->Delta[j] -
dfvl * workspace->Delta_lp_temp[j]),
CEover4 * (workspace->Delta[j] -
dfvl * workspace->Delta_lp_temp[j]),
workspace->f_ov, workspace->f_ov ); // OvCoor - 3b
Add_dDelta( system, lists, j,
CEunder4 * (1.0 - dfvl*workspace->dDelta_lp[j]) *
(bo_ij->BO_pi + bo_ij->BO_pi2),
workspace->f_un ); // UnCoor - 2a
Add_dBOpinpi2( system, lists, i, pj,
CEunder4 * (workspace->Delta[j] -
dfvl * workspace->Delta_lp_temp[j]),
CEunder4 * (workspace->Delta[j] -
dfvl * workspace->Delta_lp_temp[j]),
workspace->f_un, workspace->f_un ); // UnCoor - 2b
#endif
}
#ifdef TEST_ENERGY
//fprintf( out_control->elp, "%6d%24.15e%24.15e%24.15e\n",
//fprintf( out_control->elp, "%6d%12.4f%12.4f%12.4f\n",
// system->my_atoms[i].orig_id, workspace->nlp[i],
// e_lp, data->my_en.e_lp );
//fprintf( out_control->eov, "%6d%24.15e%24.15e\n",
fprintf( out_control->eov, "%6d%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
e_ov, data->my_en.e_ov + data->my_en.e_un );
//fprintf( out_control->eun, "%6d%24.15e%24.15e\n",
fprintf( out_control->eun, "%6d%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
e_un, data->my_en.e_ov + data->my_en.e_un );
#endif
}
}

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

@ -26,17 +26,10 @@
#include "pair_reax_c.h"
#include "reaxc_types.h"
#if defined(PURE_REAX)
#include "nonbonded.h"
#include "bond_orders.h"
#include "list.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_nonbonded.h"
#include "reaxc_bond_orders.h"
#include "reaxc_list.h"
#include "reaxc_vector.h"
#endif
void vdW_Coulomb_Energy( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
@ -56,7 +49,6 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
reax_list *far_nbrs;
// rtensor temp_rtensor, total_rtensor;
// Tallying variables:
real pe_vdw, f_tmp, delij[3];
@ -70,19 +62,17 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
e_lg = de_lg = 0.0;
for( i = 0; i < natoms; ++i ) {
if (system->my_atoms[i].type < 0) continue;
start_i = Start_Index(i, far_nbrs);
end_i = End_Index(i, far_nbrs);
orig_i = system->my_atoms[i].orig_id;
//fprintf( stderr, "i:%d, start_i: %d, end_i: %d\n", i, start_i, end_i );
for( pj = start_i; pj < end_i; ++pj ) {
nbr_pj = &(far_nbrs->select.far_nbr_list[pj]);
j = nbr_pj->nbr;
if (system->my_atoms[j].type < 0) continue;
orig_j = system->my_atoms[j].orig_id;
#if defined(PURE_REAX)
if( nbr_pj->d <= control->nonb_cut && (j < natoms || orig_i < orig_j) ) {
#elif defined(LAMMPS_REAX)
flag = 0;
if(nbr_pj->d <= control->nonb_cut) {
if (j < natoms) flag = 1;
@ -98,14 +88,11 @@ 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 ]);
/* Calculate Taper and its derivative */
// Tap = nbr_pj->Tap; -- precomputed during compte_H
Tap = workspace->Tap[7] * r_ij + workspace->Tap[6];
Tap = Tap * r_ij + workspace->Tap[5];
Tap = Tap * r_ij + workspace->Tap[4];
@ -194,18 +181,11 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
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 );
rvec_ScaledAdd( workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec );
}
else { /* NPT, iNPT or sNPT */
/* for pressure coupling, terms not related to bond order
derivatives are added directly into pressure vector/tensor */
rvec_Scale( temp, CEvd + CEclmb, nbr_pj->dvec );
rvec_ScaledAdd( workspace->f[i], -1., temp );
@ -213,46 +193,11 @@ 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] ); */
}
#ifdef TEST_ENERGY
// fprintf( out_control->evdw,
// "%12.9f%12.9f%12.9f%12.9f%12.9f%12.9f%12.9f%12.9f\n",
// workspace->Tap[7],workspace->Tap[6],workspace->Tap[5],
// workspace->Tap[4],workspace->Tap[3],workspace->Tap[2],
// workspace->Tap[1], Tap );
//fprintf( out_control->evdw, "%6d%6d%24.15e%24.15e%24.15e\n",
fprintf( out_control->evdw, "%6d%6d%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
r_ij, e_vdW, data->my_en.e_vdW );
//fprintf(out_control->ecou,"%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
fprintf( out_control->ecou, "%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
r_ij, system->my_atoms[i].q, system->my_atoms[j].q,
e_ele, data->my_en.e_ele );
#endif
#ifdef TEST_FORCES
rvec_ScaledAdd( workspace->f_vdw[i], -CEvd, nbr_pj->dvec );
rvec_ScaledAdd( workspace->f_vdw[j], +CEvd, nbr_pj->dvec );
rvec_ScaledAdd( workspace->f_ele[i], -CEclmb, nbr_pj->dvec );
rvec_ScaledAdd( workspace->f_ele[j], +CEclmb, nbr_pj->dvec );
#endif
}
}
}
#if defined(DEBUG)
fprintf( stderr, "nonbonded: ext_press (%12.6f %12.6f %12.6f)\n",
data->ext_press[0], data->ext_press[1], data->ext_press[2] );
MPI_Barrier( MPI_COMM_WORLD );
#endif
Compute_Polarization_Energy( system, data );
}
@ -284,6 +229,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
for( i = 0; i < natoms; ++i ) {
type_i = system->my_atoms[i].type;
if (type_i < 0) continue;
start_i = Start_Index(i,far_nbrs);
end_i = End_Index(i,far_nbrs);
orig_i = system->my_atoms[i].orig_id;
@ -291,11 +237,10 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
for( pj = start_i; pj < end_i; ++pj ) {
nbr_pj = &(far_nbrs->select.far_nbr_list[pj]);
j = nbr_pj->nbr;
type_j = system->my_atoms[j].type;
if (type_j < 0) continue;
orig_j = system->my_atoms[j].orig_id;
#if defined(PURE_REAX)
if( nbr_pj->d <= control->nonb_cut && (j < natoms || orig_i < orig_j) ) {
#elif defined(LAMMPS_REAX)
flag = 0;
if(nbr_pj->d <= control->nonb_cut) {
if (j < natoms) flag = 1;
@ -311,22 +256,17 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
}
if (flag) {
#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 );
t = &( LR[tmin][tmax] );
//t = &( LR[type_i][type_j] );
/* Cubic Spline Interpolation */
r = (int)(r_ij * t->inv_dx);
if( r == 0 ) ++r;
base = (real)(r+1) * t->dx;
dif = r_ij - base;
//fprintf(stderr, "r: %f, i: %d, base: %f, dif: %f\n", r, i, base, dif);
e_vdW = ((t->vdW[r].d*dif + t->vdW[r].c)*dif + t->vdW[r].b)*dif +
t->vdW[r].a;
@ -345,7 +285,6 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
t->CEclmb[r].a;
CEclmb *= system->my_atoms[i].q * system->my_atoms[j].q;
/* tally into per-atom energy */
if( system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
@ -360,8 +299,6 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
rvec_ScaledAdd( workspace->f[j], +(CEvd + CEclmb), nbr_pj->dvec );
}
else { // NPT, iNPT or sNPT
/* for pressure coupling, terms not related to bond order derivatives
are added directly into pressure vector/tensor */
rvec_Scale( temp, CEvd + CEclmb, nbr_pj->dvec );
rvec_ScaledAdd( workspace->f[i], -1., temp );
@ -370,24 +307,6 @@ void Tabulated_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 );
}
#ifdef TEST_ENERGY
//fprintf( out_control->evdw, "%6d%6d%24.15e%24.15e%24.15e\n",
fprintf( out_control->evdw, "%6d%6d%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
r_ij, e_vdW, data->my_en.e_vdW );
//fprintf(out_control->ecou,"%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
fprintf( out_control->ecou, "%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
r_ij, system->my_atoms[i].q, system->my_atoms[j].q,
e_ele, data->my_en.e_ele );
#endif
#ifdef TEST_FORCES
rvec_ScaledAdd( workspace->f_vdw[i], -CEvd, nbr_pj->dvec );
rvec_ScaledAdd( workspace->f_vdw[j], +CEvd, nbr_pj->dvec );
rvec_ScaledAdd( workspace->f_ele[i], -CEclmb, nbr_pj->dvec );
rvec_ScaledAdd( workspace->f_ele[j], +CEclmb, nbr_pj->dvec );
#endif
}
}
}
@ -404,8 +323,9 @@ void Compute_Polarization_Energy( reax_system *system, simulation_data *data )
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;
if (type_i < 0) continue;
q = system->my_atoms[i].q;
en_tmp = KCALpMOL_to_EV * (system->reax_param.sbp[type_i].chi * q +
(system->reax_param.sbp[type_i].eta / 2.) * SQR(q));
@ -507,17 +427,6 @@ void LR_vdW_Coulomb( reax_system *system, storage *workspace,
tmp = Tap / dr3gamij_3;
lr->H = EV_to_KCALpMOL * tmp;
lr->e_ele = C_ele * tmp;
// fprintf( stderr,
// "i:%d(%d), j:%d(%d), gamma:%f, Tap:%f, dr3gamij_3:%f, qi: %f, qj: %f\n",
// i, system->my_atoms[i].type, j, system->my_atoms[j].type,
// twbp->gamma, Tap, dr3gamij_3,
// system->my_atoms[i].q, system->my_atoms[j].q );
lr->CEclmb = C_ele * ( dTap - Tap * r_ij / dr3gamij_1 ) / dr3gamij_3;
// fprintf( stdout, "%d %d\t%g\t%g %g\t%g %g\t%g %g\n",
// i+1, j+1, r_ij, e_vdW, CEvd * r_ij,
// system->my_atoms[i].q, system->my_atoms[j].q, e_ele, CEclmb * r_ij );
// fprintf(stderr,"LR_Lookup: %3d %3d %5.3f-%8.5f %8.5f %8.5f %8.5f %8.5f\n",
// i, j, r_ij, lr->H, lr->e_vdW, lr->CEvd, lr->e_ele, lr->CEclmb ); */
}

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

@ -25,18 +25,10 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "reset_tools.h"
#include "list.h"
#include "tool_box.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_reset_tools.h"
#include "reaxc_list.h"
#include "reaxc_tool_box.h"
#include "reaxc_vector.h"
#endif
void Reset_Atoms( reax_system* system, control_params *control )
{
@ -47,6 +39,7 @@ void Reset_Atoms( reax_system* system, control_params *control )
if( control->hbond_cut > 0 )
for( i = 0; i < system->n; ++i ) {
atom = &(system->my_atoms[i]);
if (atom->type < 0) continue;
if( system->reax_param.sbp[ atom->type ].p_hbond == 1 )
atom->Hindex = system->numH++;
else atom->Hindex = -1;
@ -99,7 +92,6 @@ void Reset_Simulation_Data( simulation_data* data, int virial )
Reset_Energies( &data->my_en );
Reset_Energies( &data->sys_en );
Reset_Temperatures( data );
//if( virial )
Reset_Pressures( data );
}
@ -117,26 +109,6 @@ void Reset_Timing( reax_timing *rt )
rt->t_matvecs = 0;
}
#ifdef TEST_FORCES
void Reset_Test_Forces( reax_system *system, storage *workspace )
{
memset( workspace->f_ele, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_vdw, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_bo, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_be, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_lp, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_ov, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_un, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_ang, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_coa, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_pen, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_hb, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_tor, 0, system->total_cap * sizeof(rvec) );
memset( workspace->f_con, 0, system->total_cap * sizeof(rvec) );
}
#endif
void Reset_Workspace( reax_system *system, storage *workspace )
{
memset( workspace->total_bond_order, 0, system->total_cap * sizeof( real ) );
@ -144,10 +116,6 @@ void Reset_Workspace( reax_system *system, storage *workspace )
memset( workspace->CdDelta, 0, system->total_cap * sizeof( real ) );
memset( workspace->f, 0, system->total_cap * sizeof( rvec ) );
#ifdef TEST_FORCES
memset( workspace->dDelta, 0, sizeof(rvec) * system->total_cap );
Reset_Test_Forces( system, workspace );
#endif
}
@ -205,10 +173,6 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
}
}
// fprintf( stderr, "p%d: n:%d num_intrs:%d num_H:%d\n",
// system->my_rank, hbonds->n, hbonds->num_intrs, workspace->num_H );
// MPI_Barrier( comm );
/* hbonds list */
if( control->hbond_cut > 0 && system->numH > 0 ) {
hbonds = (*lists) + HBONDS;
total_hbonds = 0;
@ -234,8 +198,6 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
}
}
}
// fprintf( stderr, "p%d: cleared hbonds\n", system->my_rank );
// MPI_Barrier( comm );
}
@ -250,9 +212,4 @@ void Reset( reax_system *system, control_params *control, simulation_data *data,
Reset_Neighbor_Lists( system, control, workspace, lists, comm );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d @ step%d: reset done\n", system->my_rank, data->step );
MPI_Barrier( comm );
#endif
}

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

@ -25,16 +25,9 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "system_props.h"
#include "tool_box.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_system_props.h"
#include "reaxc_tool_box.h"
#include "reaxc_vector.h"
#endif
void Temperature_Control( control_params *control, simulation_data *data )
{
@ -255,36 +248,8 @@ void Compute_Center_of_Mass( reax_system *system, simulation_data *data,
/* Compute the rotational energy */
data->erot_cm = 0.5 * E_CONV * rvec_Dot( data->avcm, data->amcm );
#if defined(DEBUG)
fprintf( stderr, "xcm: %24.15e %24.15e %24.15e\n",
data->xcm[0], data->xcm[1], data->xcm[2] );
fprintf( stderr, "vcm: %24.15e %24.15e %24.15e\n",
data->vcm[0], data->vcm[1], data->vcm[2] );
fprintf( stderr, "amcm: %24.15e %24.15e %24.15e\n",
data->amcm[0], data->amcm[1], data->amcm[2] );
/* fprintf( stderr, "mat: %f %f %f\n %f %f %f\n %f %f %f\n",
mat[0][0], mat[0][1], mat[0][2],
mat[1][0], mat[1][1], mat[1][2],
mat[2][0], mat[2][1], mat[2][2] );
fprintf( stderr, "inv: %g %g %g\n %g %g %g\n %g %g %g\n",
inv[0][0], inv[0][1], inv[0][2],
inv[1][0], inv[1][1], inv[1][2],
inv[2][0], inv[2][1], inv[2][2] ); */
fprintf( stderr, "avcm: %24.15e %24.15e %24.15e\n",
data->avcm[0], data->avcm[1], data->avcm[2] );
#endif
}
/* IMPORTANT: This function assumes that current kinetic energy
* the system is already computed
*
* IMPORTANT: In Klein's paper, it is stated that a dU/dV term needs
* to be added when there are long-range interactions or long-range
* corrections to short-range interactions present.
* We may want to add that for more accuracy.
*/
void Compute_Pressure(reax_system* system, control_params *control,
simulation_data* data, mpi_datatypes *mpi_data)
{
@ -308,35 +273,14 @@ void Compute_Pressure(reax_system* system, control_params *control,
rvec_Multiply( tmp, p_atom->f, tx );
rvec_Add( int_press, tmp );
#if defined(DEBUG)
fprintf( stderr, "%8d%8.2f%8.2f%8.2f",
i+1, p_atom->x[0], p_atom->x[1], p_atom->x[2] );
fprintf( stderr, "%8.2f%8.2f%8.2f",
p_atom->f[0], p_atom->f[1], p_atom->f[2] );
fprintf( stderr, "%8.2f%8.2f%8.2f\n",
int_press[0], int_press[1], int_press[2] );
#endif
}
}
/* sum up internal and external pressure */
#if defined(DEBUG)
fprintf(stderr,"p%d:p_int(%10.5f %10.5f %10.5f)p_ext(%10.5f %10.5f %10.5f)\n",
system->my_rank, int_press[0], int_press[1], int_press[2],
data->my_ext_press[0], data->my_ext_press[1], data->my_ext_press[2] );
#endif
MPI_Allreduce( int_press, data->int_press,
3, MPI_DOUBLE, MPI_SUM, mpi_data->comm_mesh3D );
MPI_Allreduce( data->my_ext_press, data->ext_press,
3, MPI_DOUBLE, MPI_SUM, mpi_data->comm_mesh3D );
#if defined(DEBUG)
fprintf( stderr, "p%d: %10.5f %10.5f %10.5f\n",
system->my_rank,
data->int_press[0], data->int_press[1], data->int_press[2] );
fprintf( stderr, "p%d: %10.5f %10.5f %10.5f\n",
system->my_rank,
data->ext_press[0], data->ext_press[1], data->ext_press[2] );
#endif
/* kinetic contribution */
data->kin_press = 2.*(E_CONV*data->sys_en.e_kin) / (3.*big_box->V*P_CONV);
@ -354,63 +298,6 @@ void Compute_Pressure(reax_system* system, control_params *control,
(( data->int_press[2] + data->ext_press[2] ) /
( big_box->box_norms[0] * big_box->box_norms[1] * P_CONV ));
/* Average pressure for the whole box */
data->iso_bar.P =
( data->tot_press[0] + data->tot_press[1] + data->tot_press[2] ) / 3.;
}
/*
void Compute_Pressure_Isotropic_Klein( reax_system* system,
simulation_data* data )
{
int i;
reax_atom *p_atom;
rvec dx;
// IMPORTANT: This function assumes that current kinetic energy and
// the center of mass of the system is already computed before.
data->iso_bar.P = 2.0 * data->my_en.e_kin;
for( i = 0; i < system->N; ++i ) {
p_atom = &( system->my_atoms[i] );
rvec_ScaledSum(dx,1.0,p_atom->x,-1.0,data->xcm);
data->iso_bar.P += ( -F_CONV * rvec_Dot(p_atom->f, dx) );
}
data->iso_bar.P /= (3.0 * system->my_box.V);
// IMPORTANT: In Klein's paper, it is stated that a dU/dV term needs
// to be added when there are long-range interactions or long-range
// corrections to short-range interactions present.
// We may want to add that for more accuracy.
}
void Compute_Pressure( reax_system* system, simulation_data* data )
{
int i;
reax_atom *p_atom;
rtensor temp;
rtensor_MakeZero( data->flex_bar.P );
for( i = 0; i < system->N; ++i ) {
p_atom = &( system->my_atoms[i] );
// Distance_on_T3_Gen( data->rcm, p_atom->x, &(system->my_box), &dx );
rvec_OuterProduct( temp, p_atom->v, p_atom->v );
rtensor_ScaledAdd( data->flex_bar.P,
system->reax_param.sbp[ p_atom->type ].mass, temp );
// rvec_OuterProduct(temp,workspace->virial_forces[i],p_atom->x); //dx);
rtensor_ScaledAdd( data->flex_bar.P, -F_CONV, temp );
}
rtensor_Scale( data->flex_bar.P, 1.0 / system->my_box.V, data->flex_bar.P );
data->iso_bar.P = rtensor_Trace( data->flex_bar.P ) / 3.0;
}
*/

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

@ -25,67 +25,13 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "tool_box.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_tool_box.h"
#endif
#if defined(PURE_REAX)
/************** taken from comm_tools.c **************/
int SumScan( int n, int me, int root, MPI_Comm comm )
{
int i, my_order, wsize;;
int *nbuf = NULL;
if( me == root ) {
MPI_Comm_size( comm, &wsize );
nbuf = (int *) calloc( wsize, sizeof(int) );
MPI_Gather( &n, 1, MPI_INT, nbuf, 1, MPI_INT, root, comm );
for( i = 0; i < wsize-1; ++i ) {
nbuf[i+1] += nbuf[i];
}
MPI_Scatter( nbuf, 1, MPI_INT, &my_order, 1, MPI_INT, root, comm );
free( nbuf );
}
else {
MPI_Gather( &n, 1, MPI_INT, nbuf, 1, MPI_INT, root, comm );
MPI_Scatter( nbuf, 1, MPI_INT, &my_order, 1, MPI_INT, root, comm );
}
return my_order;
}
void SumScanB( int n, int me, int wsize, int root, MPI_Comm comm, int *nbuf )
{
int i;
MPI_Gather( &n, 1, MPI_INT, nbuf, 1, MPI_INT, root, comm );
if( me == root ) {
for( i = 0; i < wsize-1; ++i )
nbuf[i+1] += nbuf[i];
}
MPI_Bcast( nbuf, wsize, MPI_INT, root, comm );
}
#endif
/************** taken from box.c **************/
void Transform( rvec x1, simulation_box *box, char flag, rvec x2 )
{
int i, j;
real tmp;
// printf(">x1: (%lf, %lf, %lf)\n",x1[0],x1[1],x1[2]);
if (flag > 0) {
for (i=0; i < 3; i++) {
tmp = 0.0;
@ -102,7 +48,6 @@ void Transform( rvec x1, simulation_box *box, char flag, rvec x2 )
x2[i] = tmp;
}
}
// printf(">x2: (%lf, %lf, %lf)\n", x2[0], x2[1], x2[2]);
}
@ -115,8 +60,6 @@ void Transform_to_UnitBox( rvec x1, simulation_box *box, char flag, rvec x2 )
x2[2] /= box->box_norms[2];
}
/* determine whether point p is inside the box */
void Fit_to_Periodic_Box( simulation_box *box, rvec *p )
{
int i;
@ -135,126 +78,6 @@ void Fit_to_Periodic_Box( simulation_box *box, rvec *p )
}
}
#if defined(PURE_REAX)
/* determine the touch point, tp, of a box to
its neighbor denoted by the relative coordinate rl */
inline void Box_Touch_Point( simulation_box *box, ivec rl, rvec tp )
{
int d;
for( d = 0; d < 3; ++d )
if( rl[d] == -1 )
tp[d] = box->min[d];
else if( rl[d] == 0 )
tp[d] = NEG_INF - 1.;
else
tp[d] = box->max[d];
}
/* determine whether point p is inside the box */
/* assumes orthogonal box */
inline int is_Inside_Box( simulation_box *box, rvec p )
{
if( p[0] < box->min[0] || p[0] >= box->max[0] ||
p[1] < box->min[1] || p[1] >= box->max[1] ||
p[2] < box->min[2] || p[2] >= box->max[2] )
return 0;
return 1;
}
inline int iown_midpoint( simulation_box *box, rvec p1, rvec p2 )
{
rvec midp;
midp[0] = (p1[0] + p2[0]) / 2;
midp[1] = (p1[1] + p2[1]) / 2;
midp[2] = (p1[2] + p2[2]) / 2;
if( midp[0] < box->min[0] || midp[0] >= box->max[0] ||
midp[1] < box->min[1] || midp[1] >= box->max[1] ||
midp[2] < box->min[2] || midp[2] >= box->max[2] )
return 0;
return 1;
}
/**************** from grid.c ****************/
/* finds the closest point of grid cell cj to ci.
no need to consider periodic boundary conditions as in the serial case
because the box of a process is not periodic in itself */
inline void GridCell_Closest_Point( grid_cell *gci, grid_cell *gcj,
ivec ci, ivec cj, rvec cp )
{
int d;
for( d = 0; d < 3; d++ )
if( cj[d] > ci[d] )
cp[d] = gcj->min[d];
else if ( cj[d] == ci[d] )
cp[d] = NEG_INF - 1.;
else
cp[d] = gcj->max[d];
}
inline void GridCell_to_Box_Points( grid_cell *gc, ivec rl, rvec cp, rvec fp )
{
int d;
for( d = 0; d < 3; ++d )
if( rl[d] == -1 ) {
cp[d] = gc->min[d];
fp[d] = gc->max[d];
}
else if( rl[d] == 0 ) {
cp[d] = fp[d] = NEG_INF - 1.;
}
else{
cp[d] = gc->max[d];
fp[d] = gc->min[d];
}
}
inline real DistSqr_between_Special_Points( rvec sp1, rvec sp2 )
{
int i;
real d_sqr = 0;
for( i = 0; i < 3; ++i )
if( sp1[i] > NEG_INF && sp2[i] > NEG_INF )
d_sqr += SQR( sp1[i] - sp2[i] );
return d_sqr;
}
inline real DistSqr_to_Special_Point( rvec cp, rvec x )
{
int i;
real d_sqr = 0;
for( i = 0; i < 3; ++i )
if( cp[i] > NEG_INF )
d_sqr += SQR( cp[i] - x[i] );
return d_sqr;
}
inline int Relative_Coord_Encoding( ivec c )
{
return 9 * (c[0] + 1) + 3 * (c[1] + 1) + (c[2] + 1);
}
#endif
/************** from geo_tools.c *****************/
void Make_Point( real x, real y, real z, rvec* p )
{
@ -267,14 +90,6 @@ void Make_Point( real x, real y, real z, rvec* p )
int is_Valid_Serial( storage *workspace, int serial )
{
// if( workspace->map_serials[ serial ] < 0 )
// {
// fprintf( stderr, "CONECT line includes invalid pdb serial number %d.\n",
// serial );
// fprintf( stderr, "Please correct the input file.Terminating...\n" );
// MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
// }
return SUCCESS;
}
@ -303,8 +118,6 @@ void Trim_Spaces( char *element )
element[j-i] = 0; // finalize the string
}
/************ from system_props.c *************/
struct timeval tim;
real t_end;
@ -331,8 +144,6 @@ void Update_Timing_Info( real *t_start, real *timing )
*t_start = t_end;
}
/*********** from io_tools.c **************/
int Get_Atom_Type( reax_interaction *reax_param, char *s, MPI_Comm comm )
{
int i;
@ -402,8 +213,6 @@ int Tokenize( char* s, char*** tok )
return count;
}
/***************** taken from lammps ************************/
/* safe malloc */
void *smalloc( long n, const char *name, MPI_Comm comm )
{

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

@ -25,19 +25,11 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "torsion_angles.h"
#include "bond_orders.h"
#include "list.h"
#include "tool_box.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_torsion_angles.h"
#include "reaxc_bond_orders.h"
#include "reaxc_list.h"
#include "reaxc_tool_box.h"
#include "reaxc_vector.h"
#endif
#define MIN_SINE 1e-10
@ -71,12 +63,6 @@ real Calculate_Omega( rvec dvec_ij, real r_ij,
omega = atan2( unnorm_sin_omega, unnorm_cos_omega );
/* derivatives */
/* coef for adjusments to cos_theta's */
/* rla = r_ij, rlb = r_jk, rlc = r_kl, r4 = r_li;
coshd = cos_ijk, coshe = cos_jkl;
sinhd = sin_ijk, sinhe = sin_jkl; */
htra = r_ij + cos_ijk * ( r_kl * cos_jkl - r_jk );
htrb = r_jk - r_ij * cos_ijk - r_kl * cos_jkl;
htrc = r_kl + cos_jkl * ( r_ij * cos_ijk - r_jk );
@ -87,7 +73,6 @@ real Calculate_Omega( rvec dvec_ij, real r_ij,
hnhd = r_ij * r_kl * cos_ijk * sin_jkl;
hnhe = r_ij * r_kl * sin_ijk * cos_jkl;
poem = 2.0 * r_ij * r_kl * sin_ijk * sin_jkl;
if( poem < 1e-20 ) poem = 1e-20;
@ -99,27 +84,6 @@ real Calculate_Omega( rvec dvec_ij, real r_ij,
if( arg > 1.0 ) arg = 1.0;
if( arg < -1.0 ) arg = -1.0;
/* fprintf( out_control->etor,
"%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n",
htra, htrb, htrc, hthd, hthe, hnra, hnrc, hnhd, hnhe );
fprintf( out_control->etor, "%12.6f%12.6f%12.6f\n",
dvec_ij[0]/r_ij, dvec_ij[1]/r_ij, dvec_ij[2]/r_ij );
fprintf( out_control->etor, "%12.6f%12.6f%12.6f\n",
-dvec_jk[0]/r_jk, -dvec_jk[1]/r_jk, -dvec_jk[2]/r_jk );
fprintf( out_control->etor, "%12.6f%12.6f%12.6f\n",
-dvec_kl[0]/r_kl, -dvec_kl[1]/r_kl, -dvec_kl[2]/r_kl );
fprintf( out_control->etor, "%12.6f%12.6f%12.6f%12.6f\n",
r_li, dvec_li[0], dvec_li[1], dvec_li[2] );
fprintf( out_control->etor, "%12.6f%12.6f%12.6f\n",
poem, tel, arg ); */
/* fprintf( out_control->etor, "%12.6f%12.6f%12.6f\n",
-p_ijk->dcos_dk[0]/sin_ijk, -p_ijk->dcos_dk[1]/sin_ijk,
-p_ijk->dcos_dk[2]/sin_ijk );
fprintf( out_control->etor, "%12.6f%12.6f%12.6f\n",
-p_jkl->dcos_dk[0]/sin_jkl, -p_jkl->dcos_dk[1]/sin_jkl,
-p_jkl->dcos_dk[2]/sin_jkl );*/
if( sin_ijk >= 0 && sin_ijk <= MIN_SINE ) sin_ijk = MIN_SINE;
else if( sin_ijk <= 0 && sin_ijk >= -MIN_SINE ) sin_ijk = -MIN_SINE;
if( sin_jkl >= 0 && sin_jkl <= MIN_SINE ) sin_jkl = MIN_SINE;
@ -160,7 +124,7 @@ void Torsion_Angles( reax_system *system, control_params *control,
{
int i, j, k, l, pi, pj, pk, pl, pij, plk, natoms;
int type_i, type_j, type_k, type_l;
int start_j, end_j;
int start_j, end_j, start_k, end_k;
int start_pj, end_pj, start_pk, end_pk;
int num_frb_intrs = 0;
@ -198,15 +162,10 @@ void Torsion_Angles( reax_system *system, control_params *control,
real p_cot2 = system->reax_param.gp.l[27];
reax_list *bonds = (*lists) + BONDS;
reax_list *thb_intrs = (*lists) + THREE_BODIES;
// char fname[100];
// FILE *ftor;
// Virial tallying variables
real delil[3], deljl[3], delkl[3];
real eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3];
// sprintf( fname, "tor%d.out", system->my_rank );
// ftor = fopen( fname, "w" );
real eng_tmp, f_scaler, fi_tmp[3], fj_tmp[3], fk_tmp[3];
natoms = system->n;
@ -222,16 +181,12 @@ void Torsion_Angles( reax_system *system, control_params *control,
bo_jk = &( pbond_jk->bo_data );
BOA_jk = bo_jk->BO - control->thb_cut;
/* see if there are any 3-body interactions involving j&k
where j is the central atom. Otherwise there is no point in
trying to form a 4-body interaction out of this neighborhood */
if( system->my_atoms[j].orig_id < system->my_atoms[k].orig_id &&
bo_jk->BO > control->thb_cut/*0*/ && Num_Entries(pk, thb_intrs) ) {
start_k = Start_Index(k, bonds);
end_k = End_Index(k, bonds);
pj = pbond_jk->sym_index; // pj points to j on k's list
/* do the same check as above:
are there any 3-body interactions involving k&j
where k is the central atom */
if( Num_Entries(pj, thb_intrs) ) {
type_k = system->my_atoms[k].type;
Delta_k = workspace->Delta_boc[k];
@ -249,15 +204,12 @@ void Torsion_Angles( reax_system *system, control_params *control,
exp_tor34_inv = 1.0 / (1.0 + exp_tor3_DjDk + exp_tor4_DjDk);
f11_DjDk = (2.0 + exp_tor3_DjDk) * exp_tor34_inv;
/* pick i up from j-k interaction where j is the central atom */
for( pi = start_pk; pi < end_pk; ++pi ) {
p_ijk = &( thb_intrs->select.three_body_list[pi] );
pij = p_ijk->pthb; // pij is pointer to i on j's bond_list
pbond_ij = &( bonds->select.bond_list[pij] );
bo_ij = &( pbond_ij->bo_data );
if( bo_ij->BO > control->thb_cut/*0*/ ) {
i = p_ijk->thb;
type_i = system->my_atoms[i].type;
@ -277,8 +229,6 @@ void Torsion_Angles( reax_system *system, control_params *control,
exp_tor2_ij = exp( -p_tor2 * BOA_ij );
exp_cot2_ij = exp( -p_cot2 * SQR(BOA_ij -1.5) );
/* pick l up from j-k interaction where k is the central atom */
for( pl = start_pj; pl < end_pj; ++pl ) {
p_jkl = &( thb_intrs->select.three_body_list[pl] );
l = p_jkl->thb;
@ -291,7 +241,6 @@ void Torsion_Angles( reax_system *system, control_params *control,
fbp = &(system->reax_param.fbp[type_i][type_j]
[type_k][type_l].prm[0]);
if( i != l && fbh->cnt &&
bo_kl->BO > control->thb_cut/*0*/ &&
bo_ij->BO * bo_jk->BO * bo_kl->BO > control->thb_cut/*0*/ ){
@ -371,7 +320,6 @@ void Torsion_Angles( reax_system *system, control_params *control,
1.5 * fbp->V3 * (cos2omega + 2.0 * SQR(cos_omega)));
/* end of torsion energy */
/* 4-body conjugation energy */
fn12 = exp_cot2_ij * exp_cot2_jk * exp_cot2_kl;
data->my_en.e_con += e_con =
@ -512,104 +460,6 @@ void Torsion_Angles( reax_system *system, control_params *control,
if( system->pair_ptr->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",
r_ij, r_jk, r_kl, cos_ijk, cos_jkl, sin_ijk, sin_jkl );
fprintf( out_control->etor, "%12.8f\n", dfn11 ); */
/* fprintf( out_control->etor,
"%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f\n",
CEtors2, CEtors3, CEtors4, CEtors5, CEtors6,
CEtors7, CEtors8, CEtors9 ); */
/* fprintf( out_control->etor,
"%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f\n",
htra, htrb, htrc, hthd, hthe, hnra, hnrc, hnhd, hnhe ); */
/* fprintf( out_control->etor,
"%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f\n",
CEconj1, CEconj2, CEconj3, CEconj4, CEconj5, CEconj6 ); */
/* fprintf( out_control->etor, "%12.6f%12.6f%12.6f%12.6f\n",
fbp->V1, fbp->V2, fbp->V3, fbp->p_tor1 );*/
fprintf(out_control->etor,
//"%6d%6d%6d%6d%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%6d%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,system->my_atoms[l].orig_id,
RAD2DEG(omega), BOA_jk, e_tor, data->my_en.e_tor );
fprintf(out_control->econ,
//"%6d%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,system->my_atoms[l].orig_id,
RAD2DEG(omega), BOA_ij, BOA_jk, BOA_kl,
e_con, data->my_en.e_con );
#endif
#ifdef TEST_FORCES
/* Torsion Forces */
Add_dBOpinpi2( system, lists, j, pk, CEtors2, 0.0,
workspace->f_tor, workspace->f_tor );
Add_dDelta( system, lists, j, CEtors3, workspace->f_tor );
Add_dDelta( system, lists, k, CEtors3, workspace->f_tor );
Add_dBO( system, lists, j, pij, CEtors4, workspace->f_tor );
Add_dBO( system, lists, j, pk, CEtors5, workspace->f_tor );
Add_dBO( system, lists, k, plk, CEtors6, workspace->f_tor );
rvec_ScaledAdd( workspace->f_tor[i],
CEtors7, p_ijk->dcos_dk );
rvec_ScaledAdd( workspace->f_tor[j],
CEtors7, p_ijk->dcos_dj );
rvec_ScaledAdd( workspace->f_tor[k],
CEtors7, p_ijk->dcos_di );
rvec_ScaledAdd( workspace->f_tor[j],
CEtors8, p_jkl->dcos_di );
rvec_ScaledAdd( workspace->f_tor[k],
CEtors8, p_jkl->dcos_dj );
rvec_ScaledAdd( workspace->f_tor[l],
CEtors8, p_jkl->dcos_dk );
rvec_ScaledAdd( workspace->f_tor[i],
CEtors9, dcos_omega_di );
rvec_ScaledAdd( workspace->f_tor[j],
CEtors9, dcos_omega_dj );
rvec_ScaledAdd( workspace->f_tor[k],
CEtors9, dcos_omega_dk );
rvec_ScaledAdd( workspace->f_tor[l],
CEtors9, dcos_omega_dl );
/* Conjugation Forces */
Add_dBO( system, lists, j, pij, CEconj1, workspace->f_con );
Add_dBO( system, lists, j, pk, CEconj2, workspace->f_con );
Add_dBO( system, lists, k, plk, CEconj3, workspace->f_con );
rvec_ScaledAdd( workspace->f_con[i],
CEconj4, p_ijk->dcos_dk );
rvec_ScaledAdd( workspace->f_con[j],
CEconj4, p_ijk->dcos_dj );
rvec_ScaledAdd( workspace->f_con[k],
CEconj4, p_ijk->dcos_di );
rvec_ScaledAdd( workspace->f_con[j],
CEconj5, p_jkl->dcos_di );
rvec_ScaledAdd( workspace->f_con[k],
CEconj5, p_jkl->dcos_dj );
rvec_ScaledAdd( workspace->f_con[l],
CEconj5, p_jkl->dcos_dk );
rvec_ScaledAdd( workspace->f_con[i],
CEconj6, dcos_omega_di );
rvec_ScaledAdd( workspace->f_con[j],
CEconj6, dcos_omega_dj );
rvec_ScaledAdd( workspace->f_con[k],
CEconj6, dcos_omega_dk );
rvec_ScaledAdd( workspace->f_con[l],
CEconj6, dcos_omega_dl );
#endif
} // pl check ends
} // pl loop ends
} // pi check ends
@ -618,13 +468,4 @@ void Torsion_Angles( reax_system *system, control_params *control,
} // j<k && j-k neighbor check ends
} // pk loop ends
} // j loop
#if defined(DEBUG)
fprintf( stderr, "Number of torsion angles: %d\n", num_frb_intrs );
fprintf( stderr, "Torsion Energy: %g\t Conjugation Energy: %g\n",
data->my_en.e_tor, data->my_en.e_con );
fprintf( stderr, "4body: ext_press (%12.6f %12.6f %12.6f)\n",
data->ext_press[0], data->ext_press[1], data->ext_press[2] );
#endif
}

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

@ -25,56 +25,9 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "traj.h"
#include "list.h"
#include "tool_box.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_traj.h"
#include "reaxc_list.h"
#include "reaxc_tool_box.h"
#endif
#if defined(PURE_REAX)
int Set_My_Trajectory_View( MPI_File trj, int offset, MPI_Datatype etype,
MPI_Comm comm, int my_rank, int my_n, int big_n )
{
int my_disp;
int length[3];
MPI_Aint line_len;
MPI_Aint disp[3];
MPI_Datatype type[3];
MPI_Datatype view;
/* line length inferred from etype */
MPI_Type_extent( etype, &line_len );
line_len /= sizeof(char);
/* determine where to start writing into the mpi file */
my_disp = SumScan( my_n, my_rank, MASTER_NODE, comm );
my_disp -= my_n;
/* create atom_info_view */
length[0] = 1;
length[1] = my_n;
length[2] = 1;
disp[0] = 0;
disp[1] = line_len * my_disp;
disp[2] = line_len * big_n;
type[0] = MPI_LB;
type[1] = etype;
type[2] = MPI_UB;
MPI_Type_struct( 3, length, disp, type, &view );
MPI_Type_commit( &view );
MPI_File_set_view( trj, offset, etype, view, "native", MPI_INFO_NULL );
return my_disp;
}
#endif
int Reallocate_Output_Buffer( output_controls *out_control, int req_space,
MPI_Comm comm )
@ -98,31 +51,9 @@ int Reallocate_Output_Buffer( output_controls *out_control, int req_space,
void Write_Skip_Line( output_controls *out_control, mpi_datatypes *mpi_data,
int my_rank, int skip, int num_section )
{
#if defined(PURE_REAX)
MPI_Status status;
if( out_control->traj_method == MPI_TRAJ ) {
MPI_File_set_view( out_control->trj, out_control->trj_offset,
mpi_data->header_line, mpi_data->header_line,
"native", MPI_INFO_NULL );
if( my_rank == MASTER_NODE ) {
sprintf( out_control->line, INT2_LINE, "chars_to_skip_section:",
skip, num_section );
MPI_File_write( out_control->trj, out_control->line, 1,
mpi_data->header_line, &status );
}
out_control->trj_offset += HEADER_LINE_LEN;
}
else {
if( my_rank == MASTER_NODE )
fprintf( out_control->strj, INT2_LINE,
"chars_to_skip_section:", skip, num_section );
}
#elif defined(LAMMPS_REAX)
if( my_rank == MASTER_NODE )
fprintf( out_control->strj, INT2_LINE,
"chars_to_skip_section:", skip, num_section );
#endif
}
@ -191,14 +122,6 @@ int Write_Header( reax_system *system, control_params *control,
(control->restart ? "yes" : "no") );
strncat( out_control->buffer, out_control->line, HEADER_LINE_LEN+1 );
//sprintf( out_control->line, STR_LINE, "restarted_from_file:",
// (control->restart ? control->restart_from : "NA") );
//strncat( out_control->buffer, out_control->line, HEADER_LINE_LEN+1 );
//sprintf( out_control->line, STR_LINE, "kept_restart_velocities?:",
// (control->restart ? (control->random_vel ? "no":"yes"):"NA") );
//strncat( out_control->buffer, out_control->line, HEADER_LINE_LEN+1 );
sprintf( out_control->line, STR_LINE, "write_restart_files?:",
((out_control->restart_freq > 0) ? "yes" : "no") );
strncat( out_control->buffer, out_control->line, HEADER_LINE_LEN+1 );
@ -330,25 +253,8 @@ int Write_Header( reax_system *system, control_params *control,
}
/* dump out the buffer */
#if defined(PURE_REAX)
if( out_control->traj_method == MPI_TRAJ ) {
out_control->trj_offset = 0;
Set_My_Trajectory_View( out_control->trj,
out_control->trj_offset, mpi_data->header_line,
mpi_data->world, system->my_rank,
my_hdr_lines, num_hdr_lines );
MPI_File_write_all( out_control->trj, out_control->buffer,
num_hdr_lines, mpi_data->header_line, MPI_STATUS_IGNORE );
out_control->trj_offset = (num_hdr_lines) * HEADER_LINE_LEN;
}
else {
if( system->my_rank == MASTER_NODE )
fprintf( out_control->strj, "%s", out_control->buffer );
}
#elif defined(LAMMPS_REAX)
if( system->my_rank == MASTER_NODE )
fprintf( out_control->strj, "%s", out_control->buffer );
#endif
return SUCCESS;
}
@ -387,33 +293,6 @@ int Write_Init_Desc( reax_system *system, control_params *control,
out_control->line, INIT_DESC_LEN+1 );
}
#if defined(PURE_REAX)
if( out_control->traj_method == MPI_TRAJ ) {
Set_My_Trajectory_View( out_control->trj, out_control->trj_offset,
mpi_data->init_desc_line, mpi_data->world,
me, system->n, system->bigN );
MPI_File_write( out_control->trj, out_control->buffer, system->n,
mpi_data->init_desc_line, &status );
out_control->trj_offset += system->bigN * INIT_DESC_LEN;
}
else{
if( me != MASTER_NODE )
MPI_Send( out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np * INIT_DESCS + me, mpi_data->world );
else{
buffer_len = system->n * INIT_DESC_LEN;
for( i = 0; i < np; ++i )
if( i != MASTER_NODE ) {
MPI_Recv( out_control->buffer + buffer_len, buffer_req - buffer_len,
MPI_CHAR, i, np*INIT_DESCS+i, mpi_data->world, &status );
MPI_Get_count( &status, MPI_CHAR, &cnt );
buffer_len += cnt;
}
out_control->buffer[buffer_len] = 0;
fprintf( out_control->strj, "%s", out_control->buffer );
}
}
#elif defined(LAMMPS_REAX)
if( me != MASTER_NODE )
MPI_Send( out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np * INIT_DESCS + me, mpi_data->world );
@ -429,7 +308,6 @@ int Write_Init_Desc( reax_system *system, control_params *control,
out_control->buffer[buffer_len] = 0;
fprintf( out_control->strj, "%s", out_control->buffer );
}
#endif
return SUCCESS;
}
@ -464,58 +342,7 @@ int Init_Traj( reax_system *system, control_params *control,
out_control->buffer_len = 0;
out_control->buffer = NULL;
/* fprintf( stderr, "p%d: init_traj: atom_line_len = %d " \
"bond_line_len = %d, angle_line_len = %d\n" \
"max_line = %d, max_buffer_size = %d\n",
system->my_rank, out_control->atom_line_len,
out_control->bond_line_len, out_control->angle_line_len,
MAX_TRJ_LINE_LEN, MAX_TRJ_BUFFER_SIZE ); */
/* write trajectory header and atom info, if applicable */
#if defined(PURE_REAX)
if( out_control->traj_method == MPI_TRAJ ) {
/* attemp to delete the file to get rid of remnants of previous runs */
if( system->my_rank == MASTER_NODE ) {
MPI_File_delete( fname, MPI_INFO_NULL );
}
/* open a fresh trajectory file */
if( MPI_File_open( mpi_data->world, fname,
MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL,
&(out_control->trj) ) ) {
strcpy( msg, "init_traj: unable to open trajectory file" );
return FAILURE;
}
/* build the mpi structs for trajectory */
/* header_line */
MPI_Type_contiguous( HEADER_LINE_LEN, MPI_CHAR, &(mpi_data->header_line) );
MPI_Type_commit( &(mpi_data->header_line) );
/* init_desc_line */
MPI_Type_contiguous( INIT_DESC_LEN, MPI_CHAR, &(mpi_data->init_desc_line) );
MPI_Type_commit( &(mpi_data->init_desc_line) );
/* atom */
MPI_Type_contiguous( out_control->atom_line_len, MPI_CHAR,
&(mpi_data->atom_line) );
MPI_Type_commit( &(mpi_data->atom_line) );
/* bonds */
MPI_Type_contiguous( out_control->bond_line_len, MPI_CHAR,
&(mpi_data->bond_line) );
MPI_Type_commit( &(mpi_data->bond_line) );
/* angles */
MPI_Type_contiguous( out_control->angle_line_len, MPI_CHAR,
&(mpi_data->angle_line) );
MPI_Type_commit( &(mpi_data->angle_line) );
}
else if( out_control->traj_method == REG_TRAJ) {
if( system->my_rank == MASTER_NODE )
out_control->strj = fopen( fname, "w" );
}
else {
strcpy( msg, "init_traj: unknown trajectory option" );
return FAILURE;
}
#elif defined(LAMMPS_REAX)
if( out_control->traj_method == REG_TRAJ) {
if( system->my_rank == MASTER_NODE )
out_control->strj = fopen( fname, "w" );
@ -524,20 +351,8 @@ int Init_Traj( reax_system *system, control_params *control,
strcpy( msg, "init_traj: unknown trajectory option" );
return FAILURE;
}
#endif
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: initiated trajectory\n", system->my_rank );
#endif
Write_Header( system, control, out_control, mpi_data );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: header written\n", system->my_rank );
#endif
Write_Init_Desc( system, control, out_control, mpi_data );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: atom descriptions written\n", system->my_rank );
#endif
return SUCCESS;
}
@ -661,24 +476,8 @@ int Write_Frame_Header( reax_system *system, control_params *control,
}
/* dump out the buffer */
#if defined(PURE_REAX)
if( out_control->traj_method == MPI_TRAJ ) {
Set_My_Trajectory_View( out_control->trj, out_control->trj_offset,
mpi_data->header_line, mpi_data->world,
me, my_frm_hdr_lines, num_frm_hdr_lines );
MPI_File_write_all(out_control->trj, out_control->buffer, my_frm_hdr_lines,
mpi_data->header_line, MPI_STATUS_IGNORE);
out_control->trj_offset += (num_frm_hdr_lines) * HEADER_LINE_LEN;
}
else {
if( system->my_rank == MASTER_NODE )
fprintf( out_control->strj, "%s", out_control->buffer );
}
#elif defined(LAMMPS_REAX)
if( system->my_rank == MASTER_NODE )
fprintf( out_control->strj, "%s", out_control->buffer );
#endif
return SUCCESS;
}
@ -743,33 +542,6 @@ int Write_Atoms( reax_system *system, control_params *control,
strncpy( out_control->buffer + i*line_len, out_control->line, line_len+1 );
}
#if defined(PURE_REAX)
if( out_control->traj_method == MPI_TRAJ ) {
Set_My_Trajectory_View( out_control->trj, out_control->trj_offset,
mpi_data->atom_line, mpi_data->world,
me, system->n, system->bigN );
MPI_File_write( out_control->trj, out_control->buffer, system->n,
mpi_data->atom_line, &status );
out_control->trj_offset += (system->bigN) * out_control->atom_line_len;
}
else{
if( me != MASTER_NODE )
MPI_Send( out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np*ATOM_LINES+me, mpi_data->world );
else{
buffer_len = system->n * line_len;
for( i = 0; i < np; ++i )
if( i != MASTER_NODE ) {
MPI_Recv( out_control->buffer + buffer_len, buffer_req - buffer_len,
MPI_CHAR, i, np*ATOM_LINES+i, mpi_data->world, &status );
MPI_Get_count( &status, MPI_CHAR, &cnt );
buffer_len += cnt;
}
out_control->buffer[buffer_len] = 0;
fprintf( out_control->strj, "%s", out_control->buffer );
}
}
#elif defined(LAMMPS_REAX)
if( me != MASTER_NODE )
MPI_Send( out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np*ATOM_LINES+me, mpi_data->world );
@ -785,7 +557,6 @@ int Write_Atoms( reax_system *system, control_params *control,
out_control->buffer[buffer_len] = 0;
fprintf( out_control->strj, "%s", out_control->buffer );
}
#endif
return SUCCESS;
}
@ -861,34 +632,6 @@ int Write_Bonds(reax_system *system, control_params *control, reax_list *bonds,
}
}
#if defined(PURE_REAX)
if( out_control->traj_method == MPI_TRAJ ) {
Set_My_Trajectory_View( out_control->trj, out_control->trj_offset,
mpi_data->bond_line, mpi_data->world,
me, my_bonds, num_bonds );
MPI_File_write( out_control->trj, out_control->buffer, my_bonds,
mpi_data->bond_line, &status );
out_control->trj_offset += num_bonds * line_len;
}
else{
if( me != MASTER_NODE )
MPI_Send( out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np*BOND_LINES+me, mpi_data->world );
else{
buffer_len = my_bonds * line_len;
for( i = 0; i < np; ++i )
if( i != MASTER_NODE ) {
MPI_Recv( out_control->buffer + buffer_len, buffer_req - buffer_len,
MPI_CHAR, i, np*BOND_LINES+i, mpi_data->world, &status );
MPI_Get_count( &status, MPI_CHAR, &cnt );
buffer_len += cnt;
}
out_control->buffer[buffer_len] = 0;
fprintf( out_control->strj, "%s", out_control->buffer );
}
}
#elif defined(LAMMPS_REAX)
if( me != MASTER_NODE )
MPI_Send( out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np*BOND_LINES+me, mpi_data->world );
@ -904,7 +647,6 @@ int Write_Bonds(reax_system *system, control_params *control, reax_list *bonds,
out_control->buffer[buffer_len] = 0;
fprintf( out_control->strj, "%s", out_control->buffer );
}
#endif
return SUCCESS;
}
@ -985,33 +727,6 @@ int Write_Angles( reax_system *system, control_params *control,
}
}
#if defined(PURE_REAX)
if( out_control->traj_method == MPI_TRAJ ){
Set_My_Trajectory_View( out_control->trj, out_control->trj_offset,
mpi_data->angle_line, mpi_data->world,
me, my_angles, num_angles );
MPI_File_write( out_control->trj, out_control->buffer, my_angles,
mpi_data->angle_line, &status );
out_control->trj_offset += num_angles * line_len;
}
else{
if( me != MASTER_NODE )
MPI_Send( out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np*ANGLE_LINES+me, mpi_data->world );
else{
buffer_len = my_angles * line_len;
for( i = 0; i < np; ++i )
if( i != MASTER_NODE ) {
MPI_Recv( out_control->buffer + buffer_len, buffer_req - buffer_len,
MPI_CHAR, i, np*ANGLE_LINES+i, mpi_data->world, &status );
MPI_Get_count( &status, MPI_CHAR, &cnt );
buffer_len += cnt;
}
out_control->buffer[buffer_len] = 0;
fprintf( out_control->strj, "%s", out_control->buffer );
}
}
#elif defined(LAMMPS_REAX)
if( me != MASTER_NODE )
MPI_Send( out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np*ANGLE_LINES+me, mpi_data->world );
@ -1027,7 +742,6 @@ int Write_Angles( reax_system *system, control_params *control,
out_control->buffer[buffer_len] = 0;
fprintf( out_control->strj, "%s", out_control->buffer );
}
#endif
return SUCCESS;
}
@ -1037,9 +751,6 @@ int Append_Frame( reax_system *system, control_params *control,
simulation_data *data, reax_list **lists,
output_controls *out_control, mpi_datatypes *mpi_data )
{
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: appending frame %d\n", system->my_rank, data->step );
#endif
Write_Frame_Header( system, control, data, out_control, mpi_data );
if( out_control->write_atoms )
@ -1051,9 +762,6 @@ int Append_Frame( reax_system *system, control_params *control,
if( out_control->write_angles )
Write_Angles( system, control, (*lists + BONDS), (*lists + THREE_BODIES),
out_control, mpi_data );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: appended frame %d\n", system->my_rank, data->step );
#endif
return SUCCESS;
}
@ -1061,15 +769,8 @@ int Append_Frame( reax_system *system, control_params *control,
int End_Traj( int my_rank, output_controls *out_control )
{
#if defined(PURE_REAX)
if( out_control->traj_method == MPI_TRAJ )
MPI_File_close( &(out_control->trj) );
else if( my_rank == MASTER_NODE )
fclose( out_control->strj );
#elif defined(LAMMPS_REAX)
if( my_rank == MASTER_NODE )
fclose( out_control->strj );
#endif
free( out_control->buffer );
free( out_control->line );

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

@ -40,7 +40,6 @@
/************* SOME DEFS - crucial for reax_types.h *********/
//#define PURE_REAX
#define LAMMPS_REAX
//#define DEBUG
@ -48,8 +47,8 @@
//#define TEST_ENERGY
//#define TEST_FORCES
//#define CG_PERFORMANCE
#define LOG_PERFORMANCE
#define STANDARD_BOUNDARIES
//#define LOG_PERFORMANCE
//#define STANDARD_BOUNDARIES
//#define OLD_BOUNDARIES
//#define MIDPOINT_BOUNDARIES
@ -69,8 +68,6 @@ typedef real rtensor[3][3];
typedef real rvec2[2];
typedef real rvec4[4];
// import LAMMPS' definition of tagint
typedef LAMMPS_NS::tagint rc_tagint;
typedef struct {
@ -92,7 +89,6 @@ typedef struct
int type;
int num_bonds;
int num_hbonds;
//int pad; // pad to 8-byte address boundary
char name[8];
rvec x; // position
rvec v; // velocity
@ -107,20 +103,14 @@ typedef struct
int type;
int num_bonds;
int num_hbonds;
//int pad;
rvec x; // position
} boundary_atom;
typedef struct
{
//int ncells;
//int *cnt_by_gcell;
int cnt;
//int *block;
int *index;
//MPI_Datatype out_dtype;
void *out_atoms;
} mpi_out_data;
@ -147,62 +137,12 @@ typedef struct
MPI_Datatype angle_line;
MPI_Datatype angle_view;
//MPI_Request send_req1[REAX_MAX_NBRS];
//MPI_Request send_req2[REAX_MAX_NBRS];
//MPI_Status send_stat1[REAX_MAX_NBRS];
//MPI_Status send_stat2[REAX_MAX_NBRS];
//MPI_Status recv_stat1[REAX_MAX_NBRS];
//MPI_Status recv_stat2[REAX_MAX_NBRS];
mpi_out_data out_buffers[REAX_MAX_NBRS];
void *in1_buffer;
void *in2_buffer;
} mpi_datatypes;
/* Global params mapping */
/*
l[0] = p_boc1
l[1] = p_boc2
l[2] = p_coa2
l[3] = N/A
l[4] = N/A
l[5] = N/A
l[6] = p_ovun6
l[7] = N/A
l[8] = p_ovun7
l[9] = p_ovun8
l[10] = N/A
l[11] = swa
l[12] = swb
l[13] = N/A
l[14] = p_val6
l[15] = p_lp1
l[16] = p_val9
l[17] = p_val10
l[18] = N/A
l[19] = p_pen2
l[20] = p_pen3
l[21] = p_pen4
l[22] = N/A
l[23] = p_tor2
l[24] = p_tor3
l[25] = p_tor4
l[26] = N/A
l[27] = p_cot2
l[28] = p_vdW1
l[29] = v_par30
l[30] = p_coa4
l[31] = p_ovun4
l[32] = p_ovun3
l[33] = p_val8
l[34] = N/A
l[35] = N/A
l[36] = N/A
l[37] = version number
l[38] = p_coa3
*/
typedef struct
{
int n_global;
@ -815,37 +755,8 @@ typedef struct
/* force calculations */
real *CdDelta; // coefficient of dDelta
rvec *f;
#ifdef TEST_FORCES
rvec *f_ele;
rvec *f_vdw;
rvec *f_bo;
rvec *f_be;
rvec *f_lp;
rvec *f_ov;
rvec *f_un;
rvec *f_ang;
rvec *f_coa;
rvec *f_pen;
rvec *f_hb;
rvec *f_tor;
rvec *f_con;
rvec *f_tot;
rvec *dDelta; // calculated on the fly in bond_orders.c together with bo'
int *rcounts;
int *displs;
int *id_all;
rvec *f_all;
#endif
reallocate_data realloc;
//int *num_bonds;
/* hydrogen bonds */
//int num_H, Hcap;
//int *Hindex;
//int *num_hbonds;
//int *hash;
//int *rev_hash;
} storage;
@ -879,9 +790,6 @@ typedef _reax_list reax_list;
typedef struct
{
#if defined(PURE_REAX)
MPI_File trj;
#endif
FILE *strj;
int trj_offset;
int atom_line_len;
@ -916,31 +824,6 @@ typedef struct
int debug_level;
int energy_update_freq;
#ifdef TEST_ENERGY
FILE *ebond;
FILE *elp, *eov, *eun;
FILE *eval, *epen, *ecoa;
FILE *ehb;
FILE *etor, *econ;
FILE *evdw, *ecou;
#endif
#ifdef TEST_FORCES
FILE *fbo, *fdbo;
FILE *fbond;
FILE *flp, *fov, *fun;
FILE *fang, *fcoa, *fpen;
FILE *fhb;
FILE *ftor, *fcon;
FILE *fvdw, *fele;
FILE *ftot, *fcomp;
#endif
#if defined(TEST_ENERGY) || defined(TEST_FORCES)
FILE *flist; // far neighbor list
FILE *blist; // bond list
FILE *nlist; // near neighbor list
#endif
} output_controls;
@ -987,9 +870,6 @@ extern LR_lookup_table **LR;
typedef void (*evolve_function)(reax_system*, control_params*,
simulation_data*, storage*, reax_list**,
output_controls*, mpi_datatypes* );
#if defined(PURE_REAX)
evolve_function Evolve;
#endif
typedef void (*interaction_function) (reax_system*, control_params*,
simulation_data*, storage*,

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

@ -25,20 +25,11 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "valence_angles.h"
#include "bond_orders.h"
#include "list.h"
#include "vector.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_valence_angles.h"
#include "reaxc_bond_orders.h"
#include "reaxc_list.h"
#include "reaxc_vector.h"
#endif
/* calculates the theta angle between i-j-k */
void Calculate_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
real *theta, real *cos_theta )
{
@ -49,8 +40,6 @@ void Calculate_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
(*theta) = acos( *cos_theta );
}
/* calculates the derivative of the cosine of the angle between i-j-k */
void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
rvec* dcos_theta_di,
rvec* dcos_theta_dj,
@ -75,8 +64,6 @@ void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
}
/* this is a 3-body interaction in which the main role is
played by j which sits in the middle of the other two. */
void Valence_Angles( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control )
@ -104,7 +91,6 @@ void Valence_Angles( reax_system *system, control_params *control,
real r_ij, r_jk;
real BOA_ij, BOA_jk;
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];
@ -127,8 +113,8 @@ void Valence_Angles( reax_system *system, control_params *control,
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;
if (type_j < 0) continue;
start_j = Start_Index(j, bonds);
end_j = End_Index(j, bonds);
@ -145,7 +131,6 @@ void Valence_Angles( reax_system *system, control_params *control,
prod_SBO *= exp( -temp );
}
/* modifications to match Adri's code - 09/01/09 */
if( workspace->vlpex[j] >= 0 ){
vlpadj = 0;
dSBO2 = prod_SBO - 1;
@ -185,14 +170,8 @@ void Valence_Angles( reax_system *system, control_params *control,
i = pbond_ij->nbr;
r_ij = pbond_ij->d;
type_i = system->my_atoms[i].type;
// fprintf( out_control->eval, "i: %d\n", i );
/* first copy 3-body intrs from previously computed ones where i>k.
in the second for-loop below,
we compute only new 3-body intrs where i < k */
for( pk = start_j; pk < pi; ++pk ) {
// fprintf( out_control->eval, "pk: %d\n", pk );
start_pk = Start_Index( pk, thb_intrs );
end_pk = End_Index( pk, thb_intrs );
@ -213,8 +192,6 @@ void Valence_Angles( reax_system *system, control_params *control,
}
}
/* and this is the second for loop mentioned above */
for( pk = pi+1; pk < end_j; ++pk ) {
pbond_jk = &(bonds->select.bond_list[pk]);
bo_jk = &(pbond_jk->bo_data);
@ -249,23 +226,7 @@ void Valence_Angles( reax_system *system, control_params *control,
r_jk = pbond_jk->d;
thbh = &( system->reax_param.thbp[ type_i ][ type_j ][ type_k ] );
/* if( system->my_atoms[i].orig_id < system->my_atoms[k].orig_id )
fprintf( fval, "%6d %6d %6d %7.3f %7.3f %7.3f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
bo_ij->BO, bo_jk->BO, p_ijk->theta );
else
fprintf( fval, "%6d %6d %6d %7.3f %7.3f %7.3f\n",
system->my_atoms[k].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[i].orig_id,
bo_jk->BO, bo_ij->BO, p_ijk->theta ); */
for( cnt = 0; cnt < thbh->cnt; ++cnt ) {
// fprintf( out_control->eval, "%6d%6d%6d -- exists in thbp\n",
// i+1, j+1, k+1 );
if( fabs(thbh->prm[cnt].p_val1) > 0.001 ) {
thbp = &( thbh->prm[cnt] );
@ -319,7 +280,6 @@ void Valence_Angles( reax_system *system, control_params *control,
f7_ij * f7_jk * f8_Dj * expval12theta;
/* END ANGLE ENERGY*/
/* PENALTY ENERGY */
p_pen1 = thbp->p_pen1;
p_pen2 = system->reax_param.gp.l[19];
@ -346,7 +306,6 @@ void Valence_Angles( reax_system *system, control_params *control,
CEpen3 = temp * (BOA_jk - 2.0);
/* END PENALTY ENERGY */
/* COALITION ENERGY */
p_coa1 = thbp->p_coa1;
p_coa2 = system->reax_param.gp.l[2];
@ -370,7 +329,6 @@ void Valence_Angles( reax_system *system, control_params *control,
(workspace->total_bond_order[k]-BOA_jk) * e_coa;
/* END COALITION ENERGY */
/* FORCES */
bo_ij->Cdbo += (CEval1 + CEpen2 + (CEcoa1 - CEcoa4));
bo_jk->Cdbo += (CEval2 + CEpen3 + (CEcoa2 - CEcoa5));
@ -385,10 +343,6 @@ void Valence_Angles( reax_system *system, control_params *control,
temp = CUBE( temp_bo_jt );
pBOjt7 = temp * temp * temp_bo_jt;
// fprintf( out_control->eval, "%6d%12.8f\n",
// workspace->reverse_map[bonds->select.bond_list[t].nbr],
// (CEval6 * pBOjt7) );
bo_jt->Cdbo += (CEval6 * pBOjt7);
bo_jt->Cdbopi += CEval5;
bo_jt->Cdbopi2 += CEval5;
@ -400,8 +354,6 @@ void Valence_Angles( reax_system *system, control_params *control,
rvec_ScaledAdd( workspace->f[k], CEval8, p_ijk->dcos_dk );
}
else {
/* terms not related to bond order derivatives are
added directly into forces and pressure vector/tensor */
rvec_Scale( force, CEval8, p_ijk->dcos_di );
rvec_Add( workspace->f[i], force );
rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
@ -435,104 +387,6 @@ void Valence_Angles( reax_system *system, control_params *control,
if( system->pair_ptr->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",
p_val3, p_val4, BOA_ij, BOA_jk );
fprintf(out_control->eval, "%13.8f%13.8f%13.8f%13.8f%13.8f\n",
workspace->Delta_e[j], workspace->vlpex[j],
dSBO1, dSBO2, vlpadj );
fprintf( out_control->eval, "%12.8f%12.8f%12.8f%12.8f\n",
f7_ij, f7_jk, f8_Dj, expval12theta );
fprintf( out_control->eval,
"%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f%12.8f\n",
CEval1, CEval2, CEval3, CEval4,
CEval5, CEval6, CEval7, CEval8 );
fprintf( out_control->eval,
"%12.8f%12.8f%12.8f\n%12.8f%12.8f%12.8f\n%12.8f%12.8f%12.8f\n",
p_ijk->dcos_di[0]/sin_theta, p_ijk->dcos_di[1]/sin_theta,
p_ijk->dcos_di[2]/sin_theta,
p_ijk->dcos_dj[0]/sin_theta, p_ijk->dcos_dj[1]/sin_theta,
p_ijk->dcos_dj[2]/sin_theta,
p_ijk->dcos_dk[0]/sin_theta, p_ijk->dcos_dk[1]/sin_theta,
p_ijk->dcos_dk[2]/sin_theta);
fprintf( out_control->eval,
"%6d%6d%6d%15.8f%15.8f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), e_ang );*/
fprintf( out_control->eval,
//"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), theta_0, BOA_ij, BOA_jk,
e_ang, data->my_en.e_ang );
fprintf( out_control->epen,
//"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), BOA_ij, BOA_jk, e_pen,
data->my_en.e_pen );
fprintf( out_control->ecoa,
//"%6d%6d%6d%24.15e%24.15e%24.15e%24.15e%24.15e\n",
"%6d%6d%6d%12.4f%12.4f%12.4f%12.4f%12.4f\n",
system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id,
RAD2DEG(theta), BOA_ij, BOA_jk,
e_coa, data->my_en.e_coa );
#endif
#ifdef TEST_FORCES /* angle forces */
Add_dBO( system, lists, j, pi, CEval1, workspace->f_ang );
Add_dBO( system, lists, j, pk, CEval2, workspace->f_ang );
Add_dDelta( system, lists, j,
CEval3 + CEval7, workspace->f_ang );
for( t = start_j; t < end_j; ++t ) {
pbond_jt = &( bonds->select.bond_list[t] );
bo_jt = &(pbond_jt->bo_data);
temp_bo_jt = bo_jt->BO;
temp = CUBE( temp_bo_jt );
pBOjt7 = temp * temp * temp_bo_jt;
Add_dBO( system, lists, j, t, pBOjt7 * CEval6,
workspace->f_ang );
Add_dBOpinpi2( system, lists, j, t, CEval5, CEval5,
workspace->f_ang, workspace->f_ang );
}
rvec_ScaledAdd( workspace->f_ang[i], CEval8, p_ijk->dcos_di );
rvec_ScaledAdd( workspace->f_ang[j], CEval8, p_ijk->dcos_dj );
rvec_ScaledAdd( workspace->f_ang[k], CEval8, p_ijk->dcos_dk );
/* end angle forces */
/* penalty forces */
Add_dDelta( system, lists, j, CEpen1, workspace->f_pen );
Add_dBO( system, lists, j, pi, CEpen2, workspace->f_pen );
Add_dBO( system, lists, j, pk, CEpen3, workspace->f_pen );
/* end penalty forces */
/* coalition forces */
Add_dBO( system, lists, j, pi, CEcoa1 - CEcoa4,
workspace->f_coa );
Add_dBO( system, lists, j, pk, CEcoa2 - CEcoa5,
workspace->f_coa );
Add_dDelta( system, lists, j, CEcoa3, workspace->f_coa );
Add_dDelta( system, lists, i, CEcoa4, workspace->f_coa );
Add_dDelta( system, lists, k, CEcoa5, workspace->f_coa );
/* end coalition forces */
#endif
}
}
}
@ -551,16 +405,5 @@ void Valence_Angles( reax_system *system, control_params *control,
MPI_Abort( MPI_COMM_WORLD, INSUFFICIENT_MEMORY );
}
}
//fprintf( stderr,"%d: Number of angle interactions: %d\n",
// data->step, num_thb_intrs );
#if defined(DEBUG)
fprintf( stderr, "Number of angle interactions: %d\n", num_thb_intrs );
fprintf( stderr,
"Angle Energy: %g\t Penalty Energy: %g\t Coalition Energy: %g\t\n",
data->my_en.e_ang, data->my_en.e_pen, data->my_en.e_coa );
fprintf( stderr, "3body: ext_press (%12.6f %12.6f %12.6f)\n",
data->ext_press[0], data->ext_press[1], data->ext_press[2] );
#endif
}

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

@ -25,12 +25,7 @@
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
#if defined(PURE_REAX)
#include "vector.h"
#include "random.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_vector.h"
#endif
int Vector_isZero( real* v, int k )
{
@ -243,29 +238,14 @@ int rvec_isZero( rvec v )
void rvec_MakeZero( rvec v )
{
// v[0] = v[1] = v[2] = 0.0000000000000;
v[0] = v[1] = v[2] = 0.000000000000000e+00;
}
#if defined(PURE_REAX)
void rvec_Random( rvec v )
{
v[0] = Random(2.0)-1.0;
v[1] = Random(2.0)-1.0;
v[2] = Random(2.0)-1.0;
}
#endif
void rtensor_Multiply( rtensor ret, rtensor m1, rtensor m2 )
{
int i, j, k;
rtensor temp;
// check if the result matrix is the same as one of m1, m2.
// if so, we cannot modify the contents of m1 or m2, so
// we have to use a temp matrix.
if( ret == m1 || ret == m2 )
{
for( i = 0; i < 3; ++i )
@ -298,8 +278,6 @@ void rtensor_MatVec( rvec ret, rtensor m, rvec v )
int i;
rvec temp;
// if ret is the same vector as v, we cannot modify the
// contents of v until all computation is finished.
if( ret == v )
{
for( i = 0; i < 3; ++i )