compute_bond_local.cpp
d3355b03dc0elammm-master

compute_bond_local.cpp
View Options

			/* ----------------------------------------------------------------------
			LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
			http://lammps.sandia.gov, Sandia National Laboratories
			Steve Plimpton, sjplimp@sandia.gov

			Copyright (2003) Sandia Corporation. Under the terms of Contract
			DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
			certain rights in this software. This software is distributed under
			the GNU General Public License.

			See the README file in the top-level LAMMPS directory.
			------------------------------------------------------------------------- */

			#include "math.h"
			#include "string.h"
			#include "compute_bond_local.h"
			#include "atom.h"
			#include "atom_vec.h"
			#include "update.h"
			#include "domain.h"
			#include "force.h"
			#include "bond.h"
			#include "memory.h"
			#include "error.h"

			using namespace LAMMPS_NS;

			#define DELTA 10000

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

			ComputeBondLocal::ComputeBondLocal(LAMMPS lmp, int narg, char *arg) :
			Compute(lmp, narg, arg)
			{
			if (narg < 4) error->all(FLERR,"Illegal compute bond/local command");

			if (atom->avec->bonds_allow == 0)
			error->all(FLERR,"Compute bond/local used when bonds are not allowed");

			local_flag = 1;
			nvalues = narg - 3;
			if (nvalues == 1) size_local_cols = 0;
			else size_local_cols = nvalues;

			dflag = eflag = -1;
			nvalues = 0;

			int i;
			for (int iarg = 3; iarg < narg; iarg++) {
			i = iarg-3;
			if (strcmp(arg[iarg],"dist") == 0) dflag = nvalues++;
			else if (strcmp(arg[iarg],"eng") == 0) eflag = nvalues++;
			else error->all(FLERR,"Invalid keyword in compute bond/local command");
			}

			nmax = 0;
			vector = NULL;
			array = NULL;
			}

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

			ComputeBondLocal::~ComputeBondLocal()
			{
			memory->destroy(vector);
			memory->destroy(array);
			}

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

			void ComputeBondLocal::init()
			{
			if (force->bond == NULL)
			error->all(FLERR,"No bond style is defined for compute bond/local");

			// do initial memory allocation so that memory_usage() is correct

			ncount = compute_bonds(0);
			if (ncount > nmax) reallocate(ncount);
			size_local_rows = ncount;
			}

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

			void ComputeBondLocal::compute_local()
			{
			invoked_local = update->ntimestep;

			// count local entries and compute bond info

			ncount = compute_bonds(0);
			if (ncount > nmax) reallocate(ncount);
			size_local_rows = ncount;
			ncount = compute_bonds(1);
			}

			/* ----------------------------------------------------------------------
			count bonds and compute bond info on this proc
			only count bond once if newton_bond is off
			all atoms in interaction must be in group
			all atoms in interaction must be known to proc
			if bond is deleted (type = 0), do not count
			if bond is turned off (type < 0), still count
			if flag is set, compute requested info about bond
			if bond is turned off (type < 0), energy = 0.0
			------------------------------------------------------------------------- */

			int ComputeBondLocal::compute_bonds(int flag)
			{
			int i,m,n,atom1,atom2;
			double delx,dely,delz,rsq;
			double dbuf,ebuf;

			double **x = atom->x;
			int *num_bond = atom->num_bond;
			int **bond_atom = atom->bond_atom;
			int **bond_type = atom->bond_type;
			int *tag = atom->tag;
			int *mask = atom->mask;
			int nlocal = atom->nlocal;
			int newton_bond = force->newton_bond;

			if (flag) {
			if (nvalues == 1) {
			if (dflag >= 0) dbuf = vector;
			if (eflag >= 0) ebuf = vector;
			} else {
			if (dflag >= 0) dbuf = &array[0][dflag];
			if (eflag >= 0) ebuf = &array[0][eflag];
			}
			}

			Bond *bond = force->bond;

			m = n = 0;
			for (atom1 = 0; atom1 < nlocal; atom1++) {
			if (!(mask[atom1] & groupbit)) continue;
			for (i = 0; i < num_bond[atom1]; i++) {
			atom2 = atom->map(bond_atom[atom1][i]);
			if (atom2 < 0 \|\| !(mask[atom2] & groupbit)) continue;
			if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
			if (bond_type[atom1][i] == 0) continue;

			if (flag) {
			delx = x[atom1][0] - x[atom2][0];
			dely = x[atom1][1] - x[atom2][1];
			delz = x[atom1][2] - x[atom2][2];
			domain->minimum_image(delx,dely,delz);
			rsq = delxdelx + delydely + delz*delz;
			if (dflag >= 0) dbuf[n] = sqrt(rsq);
			if (eflag >= 0) {
			if (bond_type[atom1][i] > 0)
			ebuf[n] = bond->single(bond_type[atom1][i],rsq,atom1,atom2);
			else ebuf[n] = 0.0;
			}
			n += nvalues;
			}

			m++;
			}
			}

			return m;
			}

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

			void ComputeBondLocal::reallocate(int n)
			{
			// grow vector or array and indices array

			while (nmax < n) nmax += DELTA;

			if (nvalues == 1) {
			memory->destroy(vector);
			memory->create(vector,nmax,"bond/local:vector");
			vector_local = vector;
			} else {
			memory->destroy(array);
			memory->create(array,nmax,nvalues,"bond/local:array");
			array_local = array;
			}
			}

			/* ----------------------------------------------------------------------
			memory usage of local data
			------------------------------------------------------------------------- */

			double ComputeBondLocal::memory_usage()
			{
			double bytes = nmaxnvalues sizeof(double);
			return bytes;
			}

lammps/src/compute_bond_local.cppd3355b03dc0elammm-master

compute_bond_local.cppView Options

lammps/src/compute_bond_local.cpp
d3355b03dc0elammm-master

compute_bond_local.cpp
View Options