Added fix/bond/create/angle class and docu (089c4ed5) · Commits · 郑智淋 / lammps

doc/src/Commands_fix.rst

+1 −1

Original line number	Diff line number	Diff line
		@@ -42,6 +42,7 @@ OPT.
		* :doc:`bocs <fix_bocs>`
		* :doc:`bond/break <fix_bond_break>`
		* :doc:`bond/create <fix_bond_create>`
		* :doc:`bond/create/angle <fix_bond_create>`
		* :doc:`bond/react <fix_bond_react>`
		* :doc:`bond/swap <fix_bond_swap>`
		* :doc:`box/relax <fix_box_relax>`
		@@ -147,7 +148,6 @@ OPT.
		* :doc:`oneway <fix_oneway>`
		* :doc:`orient/bcc <fix_orient>`
		* :doc:`orient/fcc <fix_orient>`
		* :doc:`orient/eco <fix_orient_eco>`
		* :doc:`phonon <fix_phonon>`
		* :doc:`pimd <fix_pimd>`
		* :doc:`planeforce <fix_planeforce>`

doc/src/fix_bond_create.rst

+29 −11

Original line number	Diff line number	Diff line
		@@ -2,6 +2,8 @@

		fix bond/create command
		=======================
		fix bond/create/angle command
		=============================

		Syntax
		""""""
		@@ -17,7 +19,7 @@ Syntax
		* Rmin = 2 atoms separated by less than Rmin can bond (distance units)
		* bondtype = type of created bonds
		* zero or more keyword/value pairs may be appended to args
		* keyword = iparam or jparam or prob or atype or dtype or itype
		* keyword = iparam or jparam or prob or atype or dtype or itype or aconstrain

		.. parsed-literal::

		@@ -36,6 +38,9 @@ Syntax
		dihedraltype = type of created dihedrals
		itype value = impropertype
		impropertype = type of created impropers
		aconstrain value = amin amax
		amin = minimal angle at which new bonds can be created
		amax = maximal angle at which new bonds can be created

		Examples
		""""""""
		@@ -45,6 +50,7 @@ Examples
		fix 5 all bond/create 10 1 2 0.8 1
		fix 5 all bond/create 1 3 3 0.8 1 prob 0.5 85784 iparam 2 3
		fix 5 all bond/create 1 3 3 0.8 1 prob 0.5 85784 iparam 2 3 atype 1 dtype 2
		fix 5 all bond/create/angle 10 1 2 1.122 1 aconstrain 120 180 prob 1 4928459 iparam 2 1 jparam 2 2

		Description
		"""""""""""
		@@ -110,7 +116,16 @@ actually created. The fraction setting must be a value between 0.0
		and 1.0. A uniform random number between 0.0 and 1.0 is generated and
		the eligible bond is only created if the random number < fraction.

		Any bond that is created is assigned a bond type of bondtype
		The aconstrain keyword allows to specify a minimal and maximal angle
		amin and amax between the two prospective bonding partners and a
		third particle that is already bonded to one of the two partners.
		Such a criterion can be important, for instance when new angle
		potentials are simultaneously introduced after the formation of the
		new bond. Without a restriction on the permissible angle, and for
		stiffer angle potentials very large energies can arise and lead to
		uncontrolled behavior.

		Any bond that is created is assigned a bond type of bondtype.

		When a bond is created, data structures within LAMMPS that store bond
		topology are updated to reflect the creation. If the bond is part of
		@@ -154,13 +169,13 @@ of type angletype\ , with parameters assigned by the corresponding
		.. note::

		LAMMPS stores and maintains a data structure with a list of the
		first, second, and third neighbors of each atom (within the bond topology of
		1st, 2nd, and 3rd neighbors of each atom (within the bond topology of
		the system) for use in weighting pairwise interactions for bonded
		atoms. Note that adding a single bond always adds a new first neighbor
		but may also induce \many\ new second and third neighbors, depending on the
		atoms. Note that adding a single bond always adds a new 1st neighbor
		but may also induce \many\ new 2nd and 3rd neighbors, depending on the
		molecular topology of your system. The "extra special per atom"
		parameter must typically be set to allow for the new maximum total
		size (first + second + third neighbors) of this per-atom list. There are 2
		size (1st + 2nd + 3rd neighbors) of this per-atom list. There are 2
		ways to do this. See the :doc:`read_data <read_data>` or
		:doc:`create_box <create_box>` commands for details.

		@@ -172,12 +187,12 @@ of type angletype\ , with parameters assigned by the corresponding
		considered for pairwise interactions, using the weighting rules set by
		the :doc:`special_bonds <special_bonds>` command. Consider a new bond
		created between atoms I,J. If J has a bonded neighbor K, then K
		becomes a second neighbor of I. Even if the atype keyword is not used
		becomes a 2nd neighbor of I. Even if the atype keyword is not used
		to create angle I-J-K, the pairwise interaction between I and K will
		be potentially turned off or weighted by the 1-3 weighting specified
		by the :doc:`special_bonds <special_bonds>` command. This is the case
		even if the "angle yes" option was used with that command. The same
		is true for third neighbors (1-4 interactions), the dtype keyword, and
		is true for 3rd neighbors (1-4 interactions), the dtype keyword, and
		the "dihedral yes" option used with the
		:doc:`special_bonds <special_bonds>` command.

		@@ -237,8 +252,11 @@ Restrictions
		""""""""""""

		This fix is part of the MC package. It is only enabled if LAMMPS was
		built with that package. See the :doc:`Build package <Build_package>`
		doc page for more info.
		built with that package. See the :doc:`Build package <Build_package>` doc page for more info.

		The aconstrain keyword is only available
		when LAMMPS was built with the FixBondCreateAngle class.


		Related commands
		""""""""""""""""

src/MC/fix_bond_create.cpp

+25 −0

Original line number	Diff line number	Diff line
		@@ -79,6 +79,11 @@ FixBondCreate::FixBondCreate(LAMMPS lmp, int narg, char *arg) :
		int seed = 12345;
		atype = dtype = itype = 0;

		constrainflag = 0;
		constrainpass = 0;
		amin = 0;
		amax = 180;

		int iarg = 8;
		while (iarg < narg) {
		if (strcmp(arg[iarg],"iparam") == 0) {
		@@ -120,6 +125,20 @@ FixBondCreate::FixBondCreate(LAMMPS lmp, int narg, char *arg) :
		itype = force->inumeric(FLERR,arg[iarg+1]);
		if (itype < 0) error->all(FLERR,"Illegal fix bond/create command");
		iarg += 2;
		} else if (strcmp(arg[iarg],"aconstrain") == 0) {
		if (iarg+3 > narg) error->all(FLERR,"Illegal fix bond/create command");
		amin = force->numeric(FLERR,arg[iarg+1]);
		amax = force->inumeric(FLERR,arg[iarg+2]);
		if (amin >= amax)
		error->all(FLERR,"Illegal fix bond/create command");
		if (amin < 0 \|\| amin > 180)
		error->all(FLERR,"Illegal fix bond/create command");
		if (amax < 0 \|\| amax > 180)
		error->all(FLERR,"Illegal fix bond/create command");
		amin = (3.14159265358979323846/180.0) * amin;
		amax = (3.14159265358979323846/180.0) * amax;
		constrainflag = 1;
		iarg += 3;
		} else error->all(FLERR,"Illegal fix bond/create command");
		}

		@@ -434,6 +453,11 @@ void FixBondCreate::post_integrate()
		rsq = delxdelx + delydely + delz*delz;
		if (rsq >= cutsq) continue;

		if (constrainflag) {
		constrainpass = constrain(i,j,amin,amax);
		if (!constrainpass) continue;
		}

		if (rsq < distsq[i]) {
		partner[i] = tag[j];
		distsq[i] = rsq;
		@@ -442,6 +466,7 @@ void FixBondCreate::post_integrate()
		partner[j] = tag[i];
		distsq[j] = rsq;
		}

		}
		}

src/MC/fix_bond_create.h

+7 −2

Original line number	Diff line number	Diff line
		@@ -27,7 +27,7 @@ namespace LAMMPS_NS {
		class FixBondCreate : public Fix {
		public:
		FixBondCreate(class LAMMPS , int, char *);
		~FixBondCreate();
		virtual ~FixBondCreate();
		int setmask();
		void init();
		void init_list(int, class NeighList *);
		@@ -46,15 +46,18 @@ class FixBondCreate : public Fix {
		double compute_vector(int);
		double memory_usage();

		private:
		protected:
		int me;
		int iatomtype,jatomtype;
		int btype,seed;
		int imaxbond,jmaxbond;
		int inewtype,jnewtype;
		int constrainflag,constrainpass;
		double amin,amax;
		double cutsq,fraction;
		int atype,dtype,itype;
		int angleflag,dihedralflag,improperflag;

		int overflow;
		tagint lastcheck;

		@@ -84,6 +87,8 @@ class FixBondCreate : public Fix {
		void create_impropers(int);
		int dedup(int, int, tagint *);

		virtual int constrain(int, int, double, double) {return 1;}

		// DEBUG

		void print_bb();

src/MC/fix_bond_create_angle.cpp

0 → 100644

+145 −0

Original line number	Diff line number	Diff line
		/* ----------------------------------------------------------------------
		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.
		------------------------------------------------------------------------- */
		/* ----------------------------------------------------------------------
		Contributing authors:
		Joao Gregorio, Oliver Henrich (University of Strathclyde, Glasgow, UK)
		------------------------------------------------------------------------- */

		#include "fix_bond_create_angle.h"
		#include "atom.h"
		#include <cmath>

		using namespace LAMMPS_NS;

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

		FixBondCreateAngle::FixBondCreateAngle(LAMMPS lmp, int narg, char *arg) :
		FixBondCreate(lmp, narg, arg)
		{

		}

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

		int FixBondCreateAngle::constrain(int i, int j, double amin, double amax)
		{
		double **x = atom->x;
		tagint *tag = atom->tag;
		tagint **bond_atom = atom->bond_atom;
		int *num_bond = atom->num_bond;
		int **nspecial = atom->nspecial;
		tagint **special = atom->special;
		int *mask = atom->mask;
		int *type = atom->type;

		double v1x = 0.0;
		double v1y = 0.0;
		double v1z = 0.0;
		double v2x = 0.0;
		double v2y = 0.0;
		double v2z = 0.0;

		int flag = 0;

		// pass if both atoms have no neighbors: bond is always formed

		if (nspecial[i][0] == 0 && nspecial[j][0] == 0) flag = 1;

		// pass if i has at least one neighbor and angle constraint is met

		if (nspecial[i][0] != 0 && nspecial[j][0] == 0) {

		// calculate first vector between i and j
		// calculate second vector between i and its first neighbor

		v1x = x[i][0] - x[j][0];
		v1y = x[i][1] - x[j][1];
		v1z = x[i][2] - x[j][2];
		v2x = x[i][0] - x[atom->map(special[i][0])][0];
		v2y = x[i][1] - x[atom->map(special[i][0])][1];
		v2z = x[i][2] - x[atom->map(special[i][0])][2];

		// calculate angle between both vectors
		// set flag if the angle constraint is met

		angle1 = acos((v1xv2x + v1yv2y + v1z*v2z)/
		(sqrt(v1xv1x + v1yv1y + v1zv1z)sqrt(v2xv2x + v2yv2y + v2z*v2z)));
		if (amin <= angle1 && angle1 <= amax) flag = 1;
		}

		// pass if j has at least one neighbor and angle constraint is met

		if (nspecial[j][0] != 0 && nspecial[i][0] == 0) {

		// calculate first vector between j and i
		// calculate second vector between j and its first neighbor

		v1x = x[j][0] - x[i][0];
		v1y = x[j][1] - x[i][1];
		v1z = x[j][2] - x[i][2];
		v2x = x[j][0] - x[atom->map(special[j][0])][0];
		v2y = x[j][1] - x[atom->map(special[j][0])][1];
		v2z = x[j][2] - x[atom->map(special[j][0])][2];

		// calculate angle between both vectors
		// set flag if angle constraint is met

		angle1 = acos((v1xv2x + v1yv2y + v1z*v2z) /
		(sqrt(v1xv1x + v1yv1y + v1zv1z)sqrt(v2xv2x + v2yv2y + v2z*v2z)));
		if (amin <= angle1 && angle1 <= amax) flag = 1;
		}

		// pass if both i and j have at least one neighbor and angle constraint is met

		if (nspecial[i][0] != 0 && nspecial[j][0] != 0) {

		// Calculate 1st angle
		// calculate first vector between i and j
		// calculate second vector between i and its first neighbor

		v1x = x[i][0] - x[j][0];
		v1y = x[i][1] - x[j][1];
		v1z = x[i][2] - x[j][2];
		v2x = x[i][0] - x[atom->map(special[i][0])][0];
		v2y = x[i][1] - x[atom->map(special[i][0])][1];
		v2z = x[i][2] - x[atom->map(special[i][0])][2];

		// calculate angle between both vectors

		angle1 = acos((v1xv2x + v1yv2y + v1z*v2z) /
		(sqrt(v1xv1x + v1yv1y + v1zv1z)sqrt(v2xv2x + v2yv2y + v2z*v2z)));

		// Calculate 2nd angle
		// calculate first vector between i and j
		// calculate second vector between i and its first neighbor

		v1x = x[j][0] - x[i][0];
		v1y = x[j][1] - x[i][1];
		v1z = x[j][2] - x[i][2];
		v2x = x[j][0] - x[atom->map(special[j][0])][0];
		v2y = x[j][1] - x[atom->map(special[j][0])][1];
		v2z = x[j][2] - x[atom->map(special[j][0])][2];

		// calculate angle between both vectors

		angle2 = acos((v1xv2x + v1yv2y + v1z*v2z) /
		(sqrt(v1xv1x + v1yv1y + v1zv1z)sqrt(v2xv2x + v2yv2y + v2z*v2z)));

		// set flag if the angle constraint is met

		if ((amin <= angle1 && angle1 <= amax) && (amin <= angle2 && angle2 <= amax))
		flag = 1;
		}

		return flag;
		}

Admin message