clean up code so it can be compiled with and without OpenMP enabled regardless... (271431ab) · Commits · 郑智淋 / lammps

src/USER-OMP/fix_qeq_reax_omp.cpp

+206 −197

Original line number	Diff line number	Diff line
		@@ -18,10 +18,11 @@
		Hybrid and sub-group capabilities: Ray Shan (Sandia)
		------------------------------------------------------------------------- */

		#include "math.h"
		#include "stdio.h"
		#include "stdlib.h"
		#include "string.h"
		#include <math.h>
		#include <stdio.h>
		#include <stdlib.h>
		#include <string.h>

		#include "fix_qeq_reax_omp.h"
		#include "pair_reaxc_omp.h"
		#include "atom.h"
		@@ -321,7 +322,9 @@ void FixQEqReaxOMP::compute_H()

		// fill in the H matrix

		#if defined(_OPENMP)
		#pragma omp parallel default(shared)
		#endif
		{
		int i, j, ii, jj, mfill, jnum, flag;
		int *jlist;
		@@ -329,8 +332,9 @@ void FixQEqReaxOMP::compute_H()

		mfill = 0;

		//#pragma omp for schedule(dynamic,50)
		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for (ii = 0; ii < inum; ii++) {
		i = ilist[ii];
		if(mask[i] & groupbit) {
		@@ -390,7 +394,9 @@ void FixQEqReaxOMP::init_storage()
		if(reaxc) NN = reaxc->list->inum + reaxc->list->gnum;
		else NN = list->inum + list->gnum;

		#if defined(_OPENMP)
		#pragma omp parallel for schedule(static)
		#endif
		for (int i = 0; i < NN; i++) {
		Hdia_inv[i] = 1. / eta[atom->type[i]];
		b_s[i] = -chi[atom->type[i]];
		@@ -514,7 +520,9 @@ void FixQEqReaxOMP::init_matvec()
		if(do_aspc) {

		double m_aspc_omega = 1.0 - aspc_omega;
		#if defined(_OPENMP)
		#pragma omp parallel for schedule(dynamic,50) private(i)
		#endif
		for(int ii = 0; ii < nn; ++ii ) {
		i = ilist[ii];
		if(atom->mask[i] & groupbit) {
		@@ -540,7 +548,9 @@ void FixQEqReaxOMP::init_matvec()

		} else {

		#if defined(_OPENMP)
		#pragma omp parallel for schedule(dynamic,50) private(i)
		#endif
		for(int ii = 0; ii < nn; ++ii ) {
		i = ilist[ii];
		if(atom->mask[i] & groupbit) {
		@@ -605,7 +615,9 @@ int FixQEqReaxOMP::CG( double b, double x )
		double tmp1, tmp2;
		tmp1 = tmp2 = 0.0;

		#if defined(_OPENMP)
		#pragma omp parallel for schedule(dynamic,50) private(i) reduction(+:tmp1,tmp2)
		#endif
		for (jj = 0; jj < nn; ++jj ) {
		i = ilist[jj];
		if (atom->mask[i] & groupbit) {
		@@ -631,17 +643,23 @@ int FixQEqReaxOMP::CG( double b, double x )
		comm->reverse_comm_fix(this); //Coll_vector( q );

		tmp1 = 0.0;
		#if defined(_OPENMP)
		#pragma omp parallel
		#endif
		{

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50) private(ii) reduction(+:tmp1)
		#endif
		for (jj = 0; jj < nn; jj++) {
		ii = ilist[jj];
		if(atom->mask[ii] & groupbit) tmp1 += d[ii] * q[ii];
		}

		#if defined(_OPENMP)
		#pragma omp barrier
		#pragma omp master
		#endif
		{
		MPI_Allreduce(&tmp1, &tmp2, 1, MPI_DOUBLE, MPI_SUM, world);

		@@ -649,8 +667,10 @@ int FixQEqReaxOMP::CG( double b, double x )
		tmp1 = 0.0;
		}

		#if defined(_OPENMP)
		#pragma omp barrier
		#pragma omp for schedule(dynamic,50) private(ii) reduction(+:tmp1)
		#endif
		for (jj = 0; jj < nn; jj++) {
		ii = ilist[jj];
		if(atom->mask[ii] & groupbit) {
		@@ -671,7 +691,9 @@ int FixQEqReaxOMP::CG( double b, double x )
		sig_new = tmp2;
		beta = sig_new / sig_old;

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50) private(ii)
		#endif
		for (jj = 0; jj < nn; jj++) {
		ii = ilist[jj];
		if(atom->mask[ii] & groupbit) d[ii] = p[ii] + beta * d[ii];
		@@ -692,13 +714,19 @@ int FixQEqReaxOMP::CG( double b, double x )

		void FixQEqReaxOMP::sparse_matvec( sparse_matrix A, double x, double *b )
		{
		#if defined(_OPENMP)
		#pragma omp parallel default(shared)
		#endif
		{
		int i, j, itr_j;
		int nn, NN, ii;
		int *ilist;
		int nthreads = comm->nthreads;
		#if defined(_OPENMP)
		int tid = omp_get_thread_num();
		#else
		int tid = 0;
		#endif

		if(reaxc) {
		nn = reaxc->list->inum;
		@@ -710,26 +738,33 @@ void FixQEqReaxOMP::sparse_matvec( sparse_matrix A, double x, double *b )
		ilist = list->ilist;
		}

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (ii = 0; ii < nn; ++ii) {
		i = ilist[ii];
		if(atom->mask[i] & groupbit) b[i] = eta[ atom->type[i] ] * x[i];
		}

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (ii = nn; ii < NN; ++ii) {
		i = ilist[ii];
		if(atom->mask[i] & groupbit) b[i] = 0;
		}

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (i = 0; i < NN; ++i)
		for(int t=0; t<nthreads; t++) b_temp[t][i] = 0.0;

		// Wait for b accumulated and b_temp zeroed.
		#if defined(_OPENMP)
		#pragma omp barrier

		#pragma omp for schedule(dynamic,50)
		#endif
		for (ii = 0; ii < nn; ++ii) {
		i = ilist[ii];
		if(atom->mask[i] & groupbit) {
		@@ -743,9 +778,10 @@ void FixQEqReaxOMP::sparse_matvec( sparse_matrix A, double x, double *b )
		}

		// Wait till b_temp accumulated
		#if defined(_OPENMP)
		#pragma omp barrier

		#pragma omp for schedule(dynamic,50)
		#endif
		for (i = 0; i < NN; ++i)
		for (int t = 0; t < nthreads; ++t) b[i] += b_temp[t][i];

		@@ -772,7 +808,9 @@ void FixQEqReaxOMP::calculate_Q()

		double tmp1, tmp2;
		tmp1 = tmp2 = 0.0;
		#if defined(_OPENMP)
		#pragma omp parallel for schedule(dynamic,50) private(i) reduction(+:tmp1,tmp2)
		#endif
		for(int ii = 0; ii < nn; ii++) {
		i = ilist[ii];
		if(atom->mask[i] & groupbit) {
		@@ -792,7 +830,9 @@ void FixQEqReaxOMP::calculate_Q()

		double u = buf[0] / buf[1];

		#if defined(_OPENMP)
		#pragma omp parallel for schedule(static) private(i)
		#endif
		for (int ii = 0; ii < nn; ++ii) {
		i = ilist[ii];
		if(atom->mask[i] & groupbit) {
		@@ -814,84 +854,6 @@ void FixQEqReaxOMP::calculate_Q()

		/* ---------------------------------------------------------------------- */

		// double FixQEqReaxOMP::parallel_norm( double *v, int n )
		// {
		// int i;
		// double my_sum, norm_sqr;

		// int *ilist;

		// if (reaxc) ilist = reaxc->list->ilist;
		// else ilist = list->ilist;

		// my_sum = 0.0;
		// norm_sqr = 0.0;

		// #pragma omp parallel for schedule(static) private(i) reduction(+:my_sum)
		// for (int 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 );

		// return sqrt( norm_sqr );
		// }

		// /* ---------------------------------------------------------------------- */

		// double FixQEqReaxOMP::parallel_dot( double v1, double v2, int n )
		// {
		// int i;
		// double my_dot, res;

		// int *ilist;

		// if (reaxc) ilist = reaxc->list->ilist;
		// else ilist = list->ilist;

		// my_dot = 0.0;
		// res = 0.0;

		// #pragma omp parallel for schedule(static) private(i) reduction(+:my_dot)
		// for (int 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;
		// }

		// /* ---------------------------------------------------------------------- */

		// double FixQEqReaxOMP::parallel_vector_acc( double *v, int n )
		// {
		// int i;
		// double my_acc, res;

		// int *ilist;

		// if (reaxc) ilist = reaxc->list->ilist;
		// else ilist = list->ilist;

		// my_acc = 0.0;
		// res = 0.0;

		// #pragma omp parallel for schedule(static) private(i) reduction(+:my_acc)
		// for (int 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 );

		// return res;
		// }

		/* ---------------------------------------------------------------------- */

		void FixQEqReaxOMP::vector_sum( double* dest, double c, double* v,
		double d, double* y, int k )
		{
		@@ -901,7 +863,9 @@ void FixQEqReaxOMP::vector_sum( double* dest, double c, double* v,
		if (reaxc) ilist = reaxc->list->ilist;
		else ilist = list->ilist;

		#if defined(_OPENMP)
		#pragma omp parallel for schedule(static) private(i)
		#endif
		for (int ii=0; ii<k; ii++) {
		i = ilist[ii];
		if(atom->mask[i] & groupbit) dest[i] = c * v[i] + d * y[i];
		@@ -918,7 +882,9 @@ void FixQEqReaxOMP::vector_add( double* dest, double c, double* v, int k )
		if (reaxc) ilist = reaxc->list->ilist;
		else ilist = list->ilist;

		#if defined(_OPENMP)
		#pragma omp parallel for schedule(static) private(i)
		#endif
		for (int ii=0; ii<k; ii++) {
		i = ilist[ii];
		if(atom->mask[i] & groupbit) dest[i] += c * v[i];
		@@ -964,7 +930,9 @@ int FixQEqReaxOMP::dual_CG( double b1, double b2, double x1, double x2 )
		double tmp1, tmp2, tmp3, tmp4;
		tmp1 = tmp2 = tmp3 = tmp4 = 0.0;

		#if defined(_OPENMP)
		#pragma omp parallel for schedule(dynamic,50) private(i) reduction(+:tmp1,tmp2,tmp3,tmp4)
		#endif
		for (jj = 0; jj < nn; ++jj ) {
		i = ilist[jj];
		if (atom->mask[i] & groupbit) {
		@@ -1002,10 +970,14 @@ int FixQEqReaxOMP::dual_CG( double b1, double b2, double x1, double x2 )
		comm->reverse_comm_fix(this); //Coll_vector( q );

		tmp1 = tmp2 = 0.0;
		#if defined(_OPENMP)
		#pragma omp parallel
		#endif
		{

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50) private(ii) reduction(+:tmp1,tmp2)
		#endif
		for (jj = 0; jj < nn; jj++) {
		ii = ilist[jj];
		if(atom->mask[ii] & groupbit) {
		@@ -1015,8 +987,10 @@ int FixQEqReaxOMP::dual_CG( double b1, double b2, double x1, double x2 )
		}
		}

		#if defined(_OPENMP)
		#pragma omp barrier
		#pragma omp master
		#endif
		{
		my_buf[0] = tmp1;
		my_buf[1] = tmp2;
		@@ -1029,8 +1003,10 @@ int FixQEqReaxOMP::dual_CG( double b1, double b2, double x1, double x2 )
		tmp1 = tmp2 = 0.0;
		}

		#if defined(_OPENMP)
		#pragma omp barrier
		#pragma omp for schedule(dynamic,50) private(ii) reduction(+:tmp1,tmp2)
		#endif
		for (jj = 0; jj < nn; jj++) {
		ii = ilist[jj];
		if(atom->mask[ii] & groupbit) {
		@@ -1067,7 +1043,9 @@ int FixQEqReaxOMP::dual_CG( double b1, double b2, double x1, double x2 )
		beta_s = sig_new_s / sig_old_s;
		beta_t = sig_new_t / sig_old_t;

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50) private(ii)
		#endif
		for (jj = 0; jj < nn; jj++) {
		ii = ilist[jj];
		if(atom->mask[ii] & groupbit) {
		@@ -1131,7 +1109,9 @@ int FixQEqReaxOMP::dual_CG( double b1, double b2, double x1, double x2 )

		void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x1, double x2, double b )
		{
		#if defined(_OPENMP)
		#pragma omp parallel default(shared)
		#endif
		{
		int i, j, itr_j;
		int nn, NN, ii;
		@@ -1139,7 +1119,11 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x1, double *x2
		int indxI, indxJ;

		int nthreads = comm->nthreads;
		#if defined(_OPENMP)
		int tid = omp_get_thread_num();
		#else
		int tid = 0;
		#endif

		if (reaxc) {
		nn = reaxc->list->inum;
		@@ -1151,7 +1135,9 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x1, double *x2
		ilist = list->ilist;
		}

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (ii = 0; ii < nn; ++ii ) {
		i = ilist[ii];
		if (atom->mask[i] & groupbit) {
		@@ -1161,7 +1147,9 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x1, double *x2
		}
		}

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (ii = nn; ii < NN; ++ii ) {
		i = ilist[ii];
		if (atom->mask[i] & groupbit) {
		@@ -1171,7 +1159,9 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x1, double *x2
		}
		}

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (i = 0; i < NN; ++i) {
		indxI = 2 * i;
		for(int t=0; t<nthreads; t++) {
		@@ -1181,8 +1171,10 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x1, double *x2
		}

		// Wait for b accumulated and b_temp zeroed
		#if defined(_OPENMP)
		#pragma omp barrier
		#pragma omp for schedule(dynamic,50)
		#endif
		for (ii = 0; ii < nn; ++ii ) {
		i = ilist[ii];
		if (atom->mask[i] & groupbit) {
		@@ -1200,8 +1192,10 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x1, double *x2
		}

		// Wait till b_temp accumulated
		#if defined(_OPENMP)
		#pragma omp barrier
		#pragma omp for schedule(dynamic,50)
		#endif
		for (i = 0; i < NN; ++i) {
		indxI = 2 * i;
		for (int t = 0; t < nthreads; ++t) {
		@@ -1217,7 +1211,9 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x1, double *x2

		void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x, double *b )
		{
		#if defined(_OPENMP)
		#pragma omp parallel default(shared)
		#endif
		{
		int i, j, itr_j;
		int nn, NN, ii;
		@@ -1225,7 +1221,11 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x, double *b )
		int indxI, indxJ;

		int nthreads = comm->nthreads;
		#if defined(_OPENMP)
		int tid = omp_get_thread_num();
		#else
		int tid = 0;
		#endif

		if (reaxc) {
		nn = reaxc->list->inum;
		@@ -1237,7 +1237,9 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x, double *b )
		ilist = list->ilist;
		}

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (ii = 0; ii < nn; ++ii ) {
		i = ilist[ii];
		if (atom->mask[i] & groupbit) {
		@@ -1247,7 +1249,9 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x, double *b )
		}
		}

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (ii = nn; ii < NN; ++ii ) {
		i = ilist[ii];
		if (atom->mask[i] & groupbit) {
		@@ -1257,7 +1261,9 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x, double *b )
		}
		}

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (i = 0; i < NN; ++i) {
		indxI = 2 * i;
		for(int t=0; t<nthreads; t++) {
		@@ -1267,8 +1273,10 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x, double *b )
		}

		// Wait for b accumulated and b_temp zeroed
		#if defined(_OPENMP)
		#pragma omp barrier
		#pragma omp for schedule(dynamic,50)
		#endif
		for (ii = 0; ii < nn; ++ii ) {
		i = ilist[ii];
		if (atom->mask[i] & groupbit) {
		@@ -1286,8 +1294,10 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x, double *b )
		}

		// Wait till b_temp accumulated
		#if defined(_OPENMP)
		#pragma omp barrier
		#pragma omp for schedule(dynamic,50)
		#endif
		for (i = 0; i < NN; ++i) {
		indxI = 2 * i;
		for (int t = 0; t < nthreads; ++t) {
		@@ -1295,6 +1305,5 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix A, double x, double *b )
		b[indxI+1] += b_temp[t][indxI+1];
		}
		}

		} // omp parallel
		}

src/USER-OMP/pair_reaxc_omp.cpp

+28 −10

Original line number	Diff line number	Diff line
		@@ -49,6 +49,10 @@
		#include "reaxc_vector.h"
		#include "fix_reaxc_bonds.h"

		#if defined(_OPENMP)
		#include <omp.h>
		#endif

		using namespace LAMMPS_NS;

		#ifdef OMP_TIMING
		@@ -206,7 +210,9 @@ void PairReaxCOMP::compute(int eflag, int vflag)
		ompTimingData[COMPUTEINDEX] += (endTimeBase-startTimeBase);
		#endif

		#if defined(_OPENMP)
		#pragma omp parallel for schedule(static)
		#endif
		for(int k = 0; k < system->N; ++k) {
		num_bonds[k] = system->my_atoms[k].num_bonds;
		num_hbonds[k] = system->my_atoms[k].num_hbonds;
		@@ -269,7 +275,9 @@ void PairReaxCOMP::compute(int eflag, int vflag)
		memory->create(tmpbo,nmax,MAXSPECBOND,"pair:tmpbo");
		}

		#if defined(_OPENMP)
		#pragma omp parallel for collapse(2) schedule(static) default(shared)
		#endif
		for (i = 0; i < system->N; i ++)
		for (j = 0; j < MAXSPECBOND; j ++) {
		tmpbo[i][j] = 0.0;
		@@ -330,8 +338,11 @@ void PairReaxCOMP::init_style( )
		fix_reax = (FixReaxC *) modify->fix[modify->nfix-1];
		}

		#pragma omp parallel
		{ control->nthreads = omp_get_num_threads(); }
		#if defined(_OPENMP)
		control->nthreads = omp_get_max_threads();
		#else
		control->nthreads = 1;
		#endif
		}

		/* ---------------------------------------------------------------------- */
		@@ -418,7 +429,9 @@ void PairReaxCOMP::write_reax_atoms()
		if (system->N > system->total_cap)
		error->all(FLERR,"Too many ghost atoms");

		#if defined(_OPENMP)
		#pragma omp parallel for schedule(static) default(shared)
		#endif
		for( int i = 0; i < system->N; ++i ){
		system->my_atoms[i].orig_id = atom->tag[i];
		system->my_atoms[i].type = map[atom->type[i]];
		@@ -498,9 +511,10 @@ int PairReaxCOMP::write_reax_lists()
		num_nbrs += numneigh[i];
		}

		//#pragma omp parallel for schedule(guided) default(shared)
		#if defined(_OPENMP)
		#pragma omp parallel for schedule(dynamic,50) default(shared) \
		private(itr_i, itr_j, i, j, jlist, cutoff_sqr, num_mynbrs, d_sqr, dvec, dist)
		#endif
		for (itr_i = 0; itr_i < numall; ++itr_i) {
		i = ilist[itr_i];
		jlist = firstneigh[i];
		@@ -539,7 +553,9 @@ int PairReaxCOMP::write_reax_lists()

		void PairReaxCOMP::read_reax_forces(int vflag)
		{
		#if defined(_OPENMP)
		#pragma omp parallel for schedule(static) default(shared)
		#endif
		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];
		@@ -561,8 +577,10 @@ void PairReaxCOMP::FindBond()
		bond_data *bo_ij;
		bo_cut = 0.10;

		#if defined(_OPENMP)
		#pragma omp parallel for schedule(static) default(shared) \
		private(i, nj, pj, bo_ij, j, bo_tmp)
		#endif
		for (i = 0; i < system->n; i++) {
		nj = 0;
		for( pj = Start_Index(i, lists); pj < End_Index(i, lists); ++pj ) {

src/USER-OMP/reaxc_bond_orders_omp.cpp

+33 −6

Original line number	Diff line number	Diff line
		@@ -24,13 +24,16 @@
		<http://www.gnu.org/licenses/>.
		----------------------------------------------------------------------*/

		#include <omp.h>
		#include "pair_reaxc_omp.h"
		#include "reaxc_types.h"
		#include "reaxc_bond_orders_omp.h"
		#include "reaxc_list.h"
		#include "reaxc_vector.h"

		#if defined(_OPENMP)
		#include <omp.h>
		#endif

		using namespace LAMMPS_NS;

		/* ---------------------------------------------------------------------- */
		@@ -45,7 +48,11 @@ void Add_dBond_to_ForcesOMP( reax_system *system, int i, int pj,

		PairReaxCOMP pair_reax_ptr = static_cast<class PairReaxCOMP>(system->pair_ptr);

		#if defined(_OPENMP)
		int tid = omp_get_thread_num();
		#else
		int tid = 0;
		#endif
		ThrData *thr = pair_reax_ptr->getFixOMP()->get_thr(tid);
		long reductionOffset = (system->N * tid);

		@@ -226,7 +233,11 @@ void Add_dBond_to_Forces_NPTOMP( reax_system *system, int i, int pj, simulation_

		PairReaxCOMP pair_reax_ptr = static_cast<class PairReaxCOMP>(system->pair_ptr);

		#if defined(_OPENMP)
		int tid = omp_get_thread_num();
		#else
		int tid = 0;
		#endif
		ThrData *thr = pair_reax_ptr->getFixOMP()->get_thr(tid);
		long reductionOffset = (system->N * tid);

		@@ -428,7 +439,9 @@ void BOOMP( reax_system system, control_params control, simulation_data *data,
		int natoms = system->N;
		int nthreads = control->nthreads;

		#if defined(_OPENMP)
		#pragma omp parallel default(shared)
		#endif
		{
		int i, j, pj, type_i, type_j;
		int start_i, end_i, sym_index;
		@@ -443,10 +456,15 @@ void BOOMP( reax_system system, control_params control, simulation_data *data,
		two_body_parameters *twbp;
		bond_order_data bo_ij, bo_ji;

		#if defined(_OPENMP)
		int tid = omp_get_thread_num();

		#else
		int tid = 0;
		#endif
		/* Calculate Deltaprime, Deltaprime_boc values */
		#if defined(_OPENMP)
		#pragma omp for schedule(static)
		#endif
		for (i = 0; i < system->N; ++i) {
		type_i = system->my_atoms[i].type;
		if(type_i < 0) continue;
		@@ -459,11 +477,14 @@ void BOOMP( reax_system system, control_params control, simulation_data *data,
		}

		// Wait till initialization complete
		#if defined(_OPENMP)
		#pragma omp barrier
		#endif

		/* Corrected Bond Order calculations */
		//#pragma omp for schedule(dynamic,50)
		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for (i = 0; i < system->N; ++i) {
		type_i = system->my_atoms[i].type;
		if(type_i < 0) continue;
		@@ -625,11 +646,14 @@ void BOOMP( reax_system system, control_params control, simulation_data *data,
		}

		// Wait for bo_ij to be updated
		#if defined(_OPENMP)
		#pragma omp barrier

		#endif
		// Try to combine the following for-loop back into the for-loop above
		/-------------------------/
		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for (i = 0; i < system->N; ++i) {
		type_i = system->my_atoms[i].type;
		if(type_i < 0) continue;
		@@ -664,11 +688,14 @@ void BOOMP( reax_system system, control_params control, simulation_data *data,
		/-------------------------/

		// Need to wait for total_bond_order to be accumulated.
		#if defined(_OPENMP)
		#pragma omp barrier

		#endif
		/* Calculate some helper variables that are used at many places
		throughout force calculations */
		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for(j = 0; j < system->N; ++j ) {
		type_j = system->my_atoms[j].type;
		if(type_j < 0) continue;

src/USER-OMP/reaxc_bonds_omp.cpp

+15 −3

Original line number	Diff line number	Diff line
		@@ -25,7 +25,6 @@
		----------------------------------------------------------------------*/

		#include "pair_reaxc_omp.h"
		#include <omp.h>

		#include "reaxc_bonds_omp.h"
		#include "reaxc_bond_orders_omp.h"
		@@ -33,6 +32,10 @@
		#include "reaxc_tool_box.h"
		#include "reaxc_vector.h"

		#if defined(_OPENMP)
		#include <omp.h>
		#endif

		using namespace LAMMPS_NS;

		/* ---------------------------------------------------------------------- */
		@@ -56,7 +59,9 @@ void BondsOMP( reax_system system, control_params control,
		double gp37 = (int) system->reax_param.gp.l[37];
		double total_Ebond = 0.0;

		#if defined(_OPENMP)
		#pragma omp parallel default(shared) reduction(+: total_Ebond)
		#endif
		{
		int i, j, pj;
		int start_i, end_i;
		@@ -68,7 +73,12 @@ void BondsOMP( reax_system system, control_params control,
		single_body_parameters sbp_i, sbp_j;
		two_body_parameters *twbp;
		bond_order_data *bo_ij;

		#if defined(_OPENMP)
		int tid = omp_get_thread_num();
		#else
		int tid = 0;
		#endif
		long reductionOffset = (system->N * tid);

		class PairReaxCOMP *pair_reax_ptr;
		@@ -79,7 +89,9 @@ void BondsOMP( reax_system system, control_params control,
		system->pair_ptr->vflag_either, natoms,
		system->pair_ptr->eatom, system->pair_ptr->vatom, thr);

		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for (i = 0; i < natoms; ++i) {
		start_i = Start_Index(i, bonds);
		end_i = End_Index(i, bonds);

src/USER-OMP/reaxc_forces_omp.cpp

+50 −38

Original line number	Diff line number	Diff line
		@@ -25,7 +25,6 @@
		----------------------------------------------------------------------*/

		#include "pair_reaxc_omp.h"
		#include "omp.h"
		#include "thr_data.h"

		#include "reaxc_forces_omp.h"
		@@ -42,6 +41,10 @@
		#include "reaxc_valence_angles_omp.h"
		#include "reaxc_vector.h"

		#if defined(_OPENMP)
		#include <omp.h>
		#endif

		using namespace LAMMPS_NS;

		// Functions defined in reaxc_forces.cpp
		@@ -140,10 +143,16 @@ void Compute_Total_ForceOMP( reax_system system, control_params control,
		long totalReductionSize = system->N * nthreads;
		reax_list bonds = (lists) + BONDS;

		#if defined(_OPENMP)
		#pragma omp parallel default(shared) //default(none)
		#endif
		{
		int i, j, k, pj, pk, start_j, end_j;
		#if defined(_OPENMP)
		int tid = omp_get_thread_num();
		#else
		int tid = 0;
		#endif
		bond_order_data *bo_jk;

		class PairReaxCOMP *pair_reax_ptr;
		@@ -153,13 +162,17 @@ void Compute_Total_ForceOMP( reax_system system, control_params control,
		pair_reax_ptr->ev_setup_thr_proxy(0, 1, natoms, system->pair_ptr->eatom,
		system->pair_ptr->vatom, thr);

		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for (i = 0; i < system->N; ++i) {
		for (j = 0; j < nthreads; ++j)
		workspace->CdDelta[i] += workspace->CdDeltaReduction[system->N*j+i];
		}

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (j = 0; j < system->N; ++j) {
		start_j = Start_Index(j, bonds);
		end_j = End_Index(j, bonds);
		@@ -183,7 +196,9 @@ void Compute_Total_ForceOMP( reax_system system, control_params control,

		if(control->virial == 0) {

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (i = 0; i < system->N; ++i) {
		const int startj = Start_Index(i, bonds);
		const int endj = End_Index(i, bonds);
		@@ -194,7 +209,9 @@ void Compute_Total_ForceOMP( reax_system system, control_params control,

		} else {

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for (i = 0; i < system->N; ++i) {
		const int startj = Start_Index(i, bonds);
		const int endj = End_Index(i, bonds);
		@@ -207,14 +224,18 @@ void Compute_Total_ForceOMP( reax_system system, control_params control,

		pair_reax_ptr->reduce_thr_proxy(system->pair_ptr, 0, 1, thr);

		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for (i = 0; i < system->N; ++i) {
		for (j = 0; j < nthreads; ++j)
		rvec_Add( workspace->f[i], workspace->forceReduction[system->N*j+i] );
		}


		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for (i = 0; i < totalReductionSize; i++) {
		workspace->forceReduction[i][0] = 0;
		workspace->forceReduction[i][1] = 0;
		@@ -247,14 +268,18 @@ void Validate_ListsOMP( reax_system system, storage workspace, reax_list **lis
		reallocate_data *realloc = &(workspace->realloc);
		double saferzone = system->saferzone;

		#if defined(_OPENMP)
		#pragma omp parallel default(shared) private(i, comp, Hindex)
		#endif
		{

		/* bond list */
		if( N > 0 ) {
		bonds = *lists + BONDS;

		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for( i = 0; i < N; ++i ) {
		system->my_atoms[i].num_bonds = MAX(Num_Entries(i,bonds)*2, MIN_BONDS);

		@@ -275,7 +300,9 @@ void Validate_ListsOMP( reax_system system, storage workspace, reax_list **lis
		if( numH > 0 ) {
		hbonds = *lists + HBONDS;

		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for( i = 0; i < n; ++i ) {
		Hindex = system->my_atoms[i].Hindex;
		if( Hindex > -1 ) {
		@@ -339,18 +366,25 @@ void Init_Forces_noQEq_OMP( reax_system system, control_params control,
		/* uncorrected bond orders */
		cutoff = control->bond_cut;

		#if defined(_OPENMP)
		#pragma omp parallel default(shared) \
		private(i, atom_i, type_i, pi, start_i, end_i, sbp_i, btop_i, ibond, ihb, ihb_top, \
		j, atom_j, type_j, pj, start_j, end_j, sbp_j, nbr_pj, jbond, jhb, twbp)
		#endif
		{

		int nthreads = control->nthreads;
		#if defined(_OPENMP)
		int tid = omp_get_thread_num();
		#else
		int tid = 0;
		#endif
		long reductionOffset = system->N * tid;
		long totalReductionSize = system->N * nthreads;

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50) reduction(+ : num_bonds)
		//#pragma omp for schedule(guided) reduction(+ : num_bonds)
		#endif
		for (i = 0; i < system->N; ++i) {
		atom_i = &(system->my_atoms[i]);
		type_i = atom_i->type;
		@@ -414,7 +448,9 @@ void Init_Forces_noQEq_OMP( reax_system system, control_params control,
		int btop_j;

		// Update indices in critical section
		#if defined(_OPENMP)
		#pragma omp critical
		#endif
		{
		btop_i = End_Index( i, bonds );
		btop_j = End_Index( j, bonds );
		@@ -449,9 +485,13 @@ void Init_Forces_noQEq_OMP( reax_system system, control_params control,
		} // for(i)

		// Need to wait for all indices and tmp arrays accumulated.
		#if defined(_OPENMP)
		#pragma omp barrier
		#endif

		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50)
		#endif
		for(i=0; i<system->N; i++)
		for(int t=0; t<nthreads; t++) {
		const int indx = t*system->N + i;
		@@ -461,47 +501,22 @@ void Init_Forces_noQEq_OMP( reax_system system, control_params control,
		workspace->total_bond_order[i] += tmp_bond_order[indx];
		}

		// Not needed anymore with newest BOp_OMP()?
		//#pragma omp for schedule(guided)
		// #pragma omp for schedule(dynamic,50)
		// for (i = 0; i < system->N; ++i) {
		// start_i = Start_Index(i, bonds);
		// end_i = End_Index(i, bonds);

		// for (pi = start_i; pi < end_i; ++pi) {
		// ibond = &(bonds->select.bond_list[pi]);
		// j = ibond->nbr;

		// if (i < j) {
		// start_j = Start_Index(j, bonds);
		// end_j = End_Index(j, bonds);

		// for (pj = start_j; pj < end_j; ++pj) {
		// jbond = &(bonds->select.bond_list[pj]);

		// if (jbond->nbr == i) {
		// ibond->sym_index = pj;
		// jbond->sym_index = pi;
		// break;
		// }
		// }
		// }
		// }
		// }

		/* hydrogen bond list */
		if (control->hbond_cut > 0) {
		cutoff = control->hbond_cut;

		//#pragma omp for schedule(guided) reduction(+ : num_hbonds)
		#if defined(_OPENMP)
		#pragma omp for schedule(dynamic,50) reduction(+ : num_hbonds)
		#endif
		for (i = 0; i < system->n; ++i) {
		atom_i = &(system->my_atoms[i]);
		type_i = atom_i->type;
		sbp_i = &(system->reax_param.sbp[type_i]);
		ihb = sbp_i->p_hbond;

		#if defined(_OPENMP)
		#pragma omp critical
		#endif
		{

		if (ihb == 1 \|\| ihb == 2) {
		@@ -525,10 +540,6 @@ void Init_Forces_noQEq_OMP( reax_system system, control_params control,
		else if(j<system->n && ihb == 2 && jhb == 1) jflag = 1;

		if(iflag \|\| jflag) {

		// This critical section enforces H-bonds to be added by threads one at a time.
		// #pragma omp critical
		// {
		if(iflag) {
		ihb_top = End_Index(atom_i->Hindex, hbonds);
		Set_End_Index(atom_i->Hindex, ihb_top+1, hbonds);
		@@ -536,7 +547,6 @@ void Init_Forces_noQEq_OMP( reax_system system, control_params control,
		jhb_top = End_Index(atom_j->Hindex, hbonds);
		Set_End_Index(atom_j->Hindex, jhb_top+1, hbonds);
		}
		// } // omp critical

		if(iflag) {
		hbonds->select.hbond_list[ihb_top].nbr = j;
		@@ -561,7 +571,9 @@ void Init_Forces_noQEq_OMP( reax_system system, control_params control,
		} // if(control->hbond > 0)

		// Zero buffers for others to use as intended.
		#if defined(_OPENMP)
		#pragma omp for schedule(guided)
		#endif
		for(i=0; i<totalReductionSize; i++) {
		tmp_ddelta[i][0] = 0.0;
		tmp_ddelta[i][1] = 0.0;

Admin message