(507b038f) · Commits · 郑智淋 / lammps

src/compute_coord_atom.cpp

+4 −10

Original line number	Diff line number	Diff line
		@@ -81,11 +81,11 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS lmp, int narg, char *arg) :
		if (iorientorder < 0)
		error->all(FLERR,"Could not find compute coord/atom compute ID");
		if (strcmp(modify->compute[iorientorder]->style,"orientorder/atom") != 0)
		error->all(FLERR,"Compute coord/atom compute ID does not compute orientorder/atom");
		error->all(FLERR,"Compute coord/atom compute ID is not orientorder/atom");

		threshold = force->numeric(FLERR,arg[5]);
		if (threshold <= -1.0 \|\| threshold >= 1.0)
		error->all(FLERR,"Compute coord/atom threshold value must lie between -1 and 1");
		error->all(FLERR,"Compute coord/atom threshold not between -1 and 1");

		ncol = 1;
		typelo = new int[ncol];
		@@ -126,7 +126,7 @@ void ComputeCoordAtom::init()
		comm_forward = 2(2l+1);
		if (c_orientorder->iqlcomp < 0)
		error->all(FLERR,"Compute coord/atom requires components "
		"option in compute orientorder/atom be defined");
		"option in compute orientorder/atom");
		}

		if (force->pair == NULL)
		@@ -169,9 +169,6 @@ void ComputeCoordAtom::compute_peratom()

		invoked_peratom = update->ntimestep;

		// printf("Number of degrees %i components degree %i",nqlist,l);
		// printf("Particle \t %i \t Norm \t %g \n",0,norm[0][0]);

		// grow coordination array if necessary

		if (atom->nmax > nmax) {
		@@ -195,8 +192,6 @@ void ComputeCoordAtom::compute_peratom()
		}
		nqlist = c_orientorder->nqlist;
		int ltmp = l;
		// l = c_orientorder->qlcomp;
		if (ltmp != l) error->all(FLERR,"Debug error, ltmp != l\n");
		normv = c_orientorder->array_atom;
		comm->forward_comm_compute(this);
		}
		@@ -340,7 +335,6 @@ void ComputeCoordAtom::unpack_forward_comm(int n, int first, double *buf)
		normv[i][j] = buf[m++];
		}
		}

		}

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

src/compute_orientorder_atom.cpp

+14 −10

Original line number	Diff line number	Diff line
		@@ -73,7 +73,8 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS lmp, int narg, char *arg
		int iarg = 3;
		while (iarg < narg) {
		if (strcmp(arg[iarg],"nnn") == 0) {
		if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
		if (iarg+2 > narg)
		error->all(FLERR,"Illegal compute orientorder/atom command");
		if (strcmp(arg[iarg+1],"NULL") == 0) {
		nnn = 0;
		} else {
		@@ -91,7 +92,8 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS lmp, int narg, char *arg
		memory->destroy(qlist);
		memory->create(qlist,nqlist,"orientorder/atom:qlist");
		iarg += 2;
		if (iarg+nqlist > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
		if (iarg+nqlist > narg)
		error->all(FLERR,"Illegal compute orientorder/atom command");
		qmax = 0;
		for (int iw = 0; iw < nqlist; iw++) {
		qlist[iw] = force->numeric(FLERR,arg[iarg+iw]);
		@@ -157,11 +159,12 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
		void ComputeOrientOrderAtom::init()
		{
		if (force->pair == NULL)
		error->all(FLERR,"Compute orientorder/atom requires a pair style be defined");
		error->all(FLERR,"Compute orientorder/atom requires a "
		"pair style be defined");
		if (cutsq == 0.0) cutsq = force->pair->cutforce * force->pair->cutforce;
		else if (sqrt(cutsq) > force->pair->cutforce)
		error->all(FLERR,
		"Compute orientorder/atom cutoff is longer than pairwise cutoff");
		error->all(FLERR,"Compute orientorder/atom cutoff is "
		"longer than pairwise cutoff");

		memory->create(qnm_r,qmax,2*qmax+1,"orientorder/atom:qnm_r");
		memory->create(qnm_i,qmax,2*qmax+1,"orientorder/atom:qnm_i");
		@@ -488,7 +491,8 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
		calculate scalar distance
		------------------------------------------------------------------------- */

		double ComputeOrientOrderAtom::dist(const double r[]) {
		double ComputeOrientOrderAtom::dist(const double r[])
		{
		return sqrt(r[0]r[0] + r[1]r[1] + r[2]*r[2]);
		}

		@@ -497,8 +501,8 @@ double ComputeOrientOrderAtom::dist(const double r[]) {
		Y_l^m (theta, phi) = prefactor(l, m, cos(theta)) * exp(imphi)
		------------------------------------------------------------------------- */

		double ComputeOrientOrderAtom::
		polar_prefactor(int l, int m, double costheta) {
		double ComputeOrientOrderAtom::polar_prefactor(int l, int m, double costheta)
		{
		const int mabs = abs(m);

		double prefactor = 1.0;
		@@ -517,8 +521,8 @@ polar_prefactor(int l, int m, double costheta) {
		associated legendre polynomial
		------------------------------------------------------------------------- */

		double ComputeOrientOrderAtom::
		associated_legendre(int l, int m, double x) {
		double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x)
		{
		if (l < m) return 0.0;

		double p(1.0), pm1(0.0), pm2(0.0);

Admin message