remove redundant code from fix qeq/reax and qeq/reax/omp (0a2fe705) · Commits · 郑智淋 / lammps

src/USER-OMP/fix_qeq_reax_omp.cpp

+72 −220

Original line number	Diff line number	Diff line
		@@ -22,7 +22,6 @@
		#include <stdio.h>
		#include <stdlib.h>
		#include <string.h>

		#include "fix_qeq_reax_omp.h"
		#include "pair_reaxc_omp.h"
		#include "atom.h"
		@@ -58,36 +57,13 @@ FixQEqReaxOMP::FixQEqReaxOMP(LAMMPS lmp, int narg, char *arg) :
		{
		if (narg<8 \|\| narg>9) error->all(FLERR,"Illegal fix qeq/reax/omp command");

		nevery = force->inumeric(FLERR,arg[3]);
		swa = force->numeric(FLERR,arg[4]);
		swb = force->numeric(FLERR,arg[5]);
		tolerance = force->numeric(FLERR,arg[6]);

		// dual CG support
		dual_enabled = 0;
		if(narg == 9)
		if(strcmp(arg[8],"dual") == 0) dual_enabled = 1;
		else error->all(FLERR,"Unknown fix qeq/reax argument in location 9");

		// perform initial allocation of atom-based arrays
		// register with Atom class

		s_hist = t_hist = NULL;
		grow_arrays(atom->nmax);
		atom->add_callback(0);
		for (int i = 0; i < atom->nmax; i++ )
		for (int j = 0; j < nprev; ++j )
		s_hist[i][j] = t_hist[i][j] = 0;

		reaxc = NULL;
		reaxc = (PairReaxCOMP *) force->pair_match("reax/c/omp",1);

		b_temp = NULL;

		// ASPC: Kolafa, J. Comp. Chem., 25(3), 335 (2003)
		do_aspc = 0;
		aspc_order = 1;
		aspc_order_max = nprev - 2; // Must be consistent with nprev to store history: nprev = aspc_order + 2
		// Must be consistent with nprev to store history: nprev = aspc_order + 2
		aspc_order_max = nprev - 2;
		aspc_omega = 0.0;
		aspc_b = NULL;
		}
		@@ -106,56 +82,6 @@ void FixQEqReaxOMP::post_constructor()

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

		void FixQEqReaxOMP::pertype_parameters(char *arg)
		{
		if (strcmp(arg,"reax/c") == 0) {
		reaxflag = 1;
		Pair *pair = force->pair_match("reax/c",0);
		if (pair == NULL) error->all(FLERR,"No pair reax/c for fix qeq/reax/omp");
		int tmp;
		chi = (double *) pair->extract("chi",tmp);
		eta = (double *) pair->extract("eta",tmp);
		gamma = (double *) pair->extract("gamma",tmp);
		if (chi == NULL \|\| eta == NULL \|\| gamma == NULL)
		error->all(FLERR,
		"Fix qeq/reax/omp could not extract params from pair reax/c");
		return;
		}

		int i,itype,ntypes;
		double v1,v2,v3;
		FILE *pf;

		reaxflag = 0;
		ntypes = atom->ntypes;

		memory->create(chi,ntypes+1,"qeq/reax:chi");
		memory->create(eta,ntypes+1,"qeq/reax:eta");
		memory->create(gamma,ntypes+1,"qeq/reax:gamma");

		if (comm->me == 0) {
		if ((pf = fopen(arg,"r")) == NULL)
		error->one(FLERR,"Fix qeq/reax/omp parameter file could not be found");

		for (i = 1; i <= ntypes && !feof(pf); i++) {
		fscanf(pf,"%d %lg %lg %lg",&itype,&v1,&v2,&v3);
		if (itype < 1 \|\| itype > ntypes)
		error->one(FLERR,"Fix qeq/reax/omp invalid atom type in param file");
		chi[itype] = v1;
		eta[itype] = v2;
		gamma[itype] = v3;
		}
		if (i <= ntypes) error->one(FLERR,"Invalid param file for fix qeq/reax/omp");
		fclose(pf);
		}

		MPI_Bcast(&chi[1],ntypes,MPI_DOUBLE,0,world);
		MPI_Bcast(&eta[1],ntypes,MPI_DOUBLE,0,world);
		MPI_Bcast(&gamma[1],ntypes,MPI_DOUBLE,0,world);
		}

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

		void FixQEqReaxOMP::allocate_storage()
		{
		FixQEqReax::allocate_storage();
		@@ -177,78 +103,9 @@ void FixQEqReaxOMP::deallocate_storage()

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

		void FixQEqReaxOMP::allocate_matrix()
		{
		int i,ii,inum,m;
		int ilist, numneigh;

		int mincap;
		double safezone;

		if( reaxflag ) {
		mincap = reaxc->system->mincap;
		safezone = reaxc->system->safezone;
		} else {
		mincap = MIN_CAP;
		safezone = SAFE_ZONE;
		}

		n = atom->nlocal;
		n_cap = MAX( (int)(n * safezone), mincap );

		// determine the total space for the H matrix

		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 );

		H.n = n_cap;
		H.m = m_cap;
		memory->create(H.firstnbr,n_cap,"qeq:H.firstnbr");
		memory->create(H.numnbrs,n_cap,"qeq:H.numnbrs");
		memory->create(H.jlist,m_cap,"qeq:H.jlist");
		memory->create(H.val,m_cap,"qeq:H.val");
		}

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

		void FixQEqReaxOMP::init()
		{
		if (!atom->q_flag) error->all(FLERR,"Fix qeq/reax/omp requires atom attribute q");

		ngroup = group->count(igroup);
		if (ngroup == 0) error->all(FLERR,"Fix qeq/reax group has no atoms");

		if (!force->pair_match("reax/c/omp",1))
		error->all(FLERR,"Must use pair_style reax/c/omp with fix qeq/reax/omp");

		// need a half neighbor list w/ Newton off and ghost neighbors
		// built whenever re-neighboring occurs

		int irequest = neighbor->request(this,instance_me);
		neighbor->requests[irequest]->pair = 0;
		neighbor->requests[irequest]->fix = 1;
		neighbor->requests[irequest]->newton = 2;
		neighbor->requests[irequest]->ghost = 1;

		init_shielding();
		init_taper();

		if (strstr(update->integrate_style,"respa"))
		nlevels_respa = ((Respa *) update->integrate)->nlevels;
		FixQEqReax::init();

		// APSC setup
		if (do_aspc) {
		@@ -276,10 +133,6 @@ void FixQEqReaxOMP::init()
		n -= 1.0;
		d += 1.0;
		}

		// fprintf(stdout,"aspc_omega= %f\n",aspc_omega);
		// for(int i=0; i<aspc_order_max+2; i++) fprintf(stdout,"i= %i B= %f\n",i,aspc_b[i]);
		// error->all(FLERR,"Early Termination");
		}
		}

		@@ -444,9 +297,9 @@ void FixQEqReaxOMP::pre_force(int vflag)
		startTimeBase = endTimeBase;
		#endif

		if(dual_enabled) matvecs = dual_CG(b_s, b_t, s, t); // OMP_TIMING inside dual_CG
		else {

		if (dual_enabled) {
		matvecs = dual_CG(b_s, b_t, s, t); // OMP_TIMING inside dual_CG
		} else {
		matvecs_s = CG(b_s, s); // CG on s - parallel

		#ifdef OMP_TIMING
		@@ -469,8 +322,6 @@ void FixQEqReaxOMP::pre_force(int vflag)

		} // if (dual_enabled)

		// if(comm->me == 0) fprintf(stdout,"matvecs= %i %i\n",matvecs_s,matvecs_t);

		#ifdef OMP_TIMING
		startTimeBase = MPI_Wtime();
		#endif
		@@ -590,8 +441,8 @@ void FixQEqReaxOMP::init_matvec()

		int FixQEqReaxOMP::CG( double b, double x)
		{
		int i, ii, j, imax;
		double tmp, alpha, beta, b_norm;
		int i, ii, imax;
		double alpha, beta, b_norm;
		double sig_old, sig_new;

		double my_buf[2], buf[2];
		@@ -905,8 +756,8 @@ int FixQEqReaxOMP::dual_CG( double b1, double b2, double x1, double x2 )
		startTimeBase = MPI_Wtime();
		#endif

		int i, j, imax;
		double tmp, alpha_s, alpha_t, beta_s, beta_t, b_norm_s, b_norm_t;
		int i, imax;
		double alpha_s, alpha_t, beta_s, beta_t, b_norm_s, b_norm_t;
		double sig_old_s, sig_old_t, sig_new_s, sig_new_t;

		double my_buf[4], buf[4];
		@@ -1038,7 +889,8 @@ int FixQEqReaxOMP::dual_CG( double b1, double b2, double x1, double x2 )
		sig_new_s = buf[0];
		sig_new_t = buf[1];

		if( sqrt(sig_new_s)/b_norm_s <= tolerance \|\| sqrt(sig_new_t)/b_norm_t <= tolerance) break;
		if (sqrt(sig_new_s)/b_norm_s <= tolerance
		\|\| sqrt(sig_new_t)/b_norm_t <= tolerance) break;

		beta_s = sig_new_s / sig_old_s;
		beta_t = sig_new_t / sig_old_t;

src/USER-OMP/fix_qeq_reax_omp.h

+0 −8

Original line number	Diff line number	Diff line
		@@ -47,17 +47,13 @@ class FixQEqReaxOMP : public FixQEqReax {
		protected:
		double **b_temp;

		class PairReaxCOMP *reaxc;

		int do_aspc;
		int aspc_order, aspc_order_max;
		double aspc_omega;
		double * aspc_b;

		virtual void pertype_parameters(char*);
		virtual void allocate_storage();
		virtual void deallocate_storage();
		virtual void allocate_matrix();
		virtual void init_matvec();
		virtual void compute_H();

		@@ -65,10 +61,6 @@ class FixQEqReaxOMP : public FixQEqReax {
		virtual void sparse_matvec(sparse_matrix,double,double*);
		virtual void calculate_Q();

		/* virtual double parallel_norm( double, int ); /
		/* virtual double parallel_dot( double, double, int ); */
		/* virtual double parallel_vector_acc( double, int ); /

		virtual void vector_sum(double,double,double,double,double*,int);
		virtual void vector_add(double, double, double,int);

src/USER-REAXC/fix_qeq_reax.cpp

+89 −96

Original line number	Diff line number	Diff line
		@@ -81,8 +81,12 @@ FixQEqReax::FixQEqReax(LAMMPS lmp, int narg, char *arg) :
		strcpy(pertype_option,arg[7]);

		// dual CG support only available for USER-OMP variant
		// check for compatibility is in Fix::post_constructor()
		dual_enabled = 0;

		if (narg == 9) {
		if (strcmp(arg[8],"dual") == 0) dual_enabled = 1;
		else error->all(FLERR,"Illegal fix qeq/reax command");
		}
		shld = NULL;

		n = n_cap = 0;
		@@ -111,13 +115,16 @@ FixQEqReax::FixQEqReax(LAMMPS lmp, int narg, char *arg) :
		H.jlist = NULL;
		H.val = NULL;

		comm_forward = comm_reverse = 1;
		// dual CG support
		// Update comm sizes for this fix
		if (dual_enabled) comm_forward = comm_reverse = 2;
		else comm_forward = comm_reverse = 1;

		// perform initial allocation of atom-based arrays
		// register with Atom class

		reaxc = NULL;
		reaxc = (PairReaxC *) force->pair_match("reax/c",1);
		reaxc = (PairReaxC *) force->pair_match("reax/c",0);

		if (reaxc) {
		s_hist = t_hist = NULL;
		@@ -127,10 +134,6 @@ FixQEqReax::FixQEqReax(LAMMPS lmp, int narg, char *arg) :
		for (int j = 0; j < nprev; ++j)
		s_hist[i][j] = t_hist[i][j] = 0;
		}

		// dual CG support
		// Update comm sizes for this fix
		if (dual_enabled) comm_forward = comm_reverse = 2;
		}

		/* ---------------------------------------------------------------------- */
		@@ -165,6 +168,8 @@ FixQEqReax::~FixQEqReax()
		void FixQEqReax::post_constructor()
		{
		pertype_parameters(pertype_option);
		if (dual_enabled)
		error->all(FLERR,"Dual keyword only supported with fix qeq/reax/omp");
		}

		/* ---------------------------------------------------------------------- */
		@@ -184,11 +189,9 @@ void FixQEqReax::pertype_parameters(char *arg)
		{
		if (strcmp(arg,"reax/c") == 0) {
		reaxflag = 1;
		Pair *pair = force->pair_match("reax/c",1);
		if (pair == NULL)
		pair = force->pair_match("reax/c/kk",1);

		Pair *pair = force->pair_match("reax/c",0);
		if (pair == NULL) error->all(FLERR,"No pair reax/c for fix qeq/reax");

		int tmp;
		chi = (double *) pair->extract("chi",tmp);
		eta = (double *) pair->extract("eta",tmp);
		@@ -199,10 +202,6 @@ void FixQEqReax::pertype_parameters(char *arg)
		return;
		}

		// OMP style will use it's own pertype_parameters()
		Pair * pair = force->pair_match("reax/c/omp",1);
		if (pair) return;

		int i,itype,ntypes;
		double v1,v2,v3;
		FILE *pf;
		@@ -358,18 +357,12 @@ void FixQEqReax::reallocate_matrix()

		void FixQEqReax::init()
		{
		if (!atom->q_flag) error->all(FLERR,"Fix qeq/reax requires atom attribute q");
		if (!atom->q_flag)
		error->all(FLERR,"Fix qeq/reax requires atom attribute q");

		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
		// built whenever re-neighboring occurs

		@@ -590,12 +583,12 @@ void FixQEqReax::compute_H()
		{
		int inum, jnum, ilist, jlist, numneigh, *firstneigh;
		int i, j, ii, jj, flag;
		double **x, SMALL = 0.0001;
		double dx, dy, dz, r_sqr;
		const double SMALL = 0.0001;

		int *type = atom->type;
		tagint *tag = atom->tag;
		x = atom->x;
		double **x = atom->x;
		int *mask = atom->mask;

		if (reaxc) {

Admin message