Merge branch 'master' into write-bonus-data

   * :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>`

* :doc:`bocs <fix_bocs>` - NPT style time integration with pressure correction

* :doc:`bond/break <fix_bond_break>` - break bonds on the fly

* :doc:`bond/create <fix_bond_create>` - create bonds on the fly

* :doc:`bond/create/angle <fix_bond_create>` - create bonds on the fly with angle constraints

* :doc:`bond/react <fix_bond_react>` - apply topology changes to model reactions

* :doc:`bond/swap <fix_bond_swap>` - Monte Carlo bond swapping

* :doc:`box/relax <fix_box_relax>` - relax box size during energy minimization

fix bond/create command

=======================

fix bond/create/angle command

=============================

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::

         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

""""""""

   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

"""""""""""

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 is only available with the fix

bond/create/angle command.  It 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 when new angles are defined together

with the formation of a 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

**Restart, fix_modify, output, run start/stop, minimize info:**

No information about this fix is written to :doc:`binary restart files <restart>`.  None of the :doc:`fix_modify <fix_modify>` options

are relevant to this fix.

No information about this fix is written to :doc:`binary restart files

<restart>`.  None of the :doc:`fix_modify <fix_modify>` options are

relevant to this fix.

This fix computes two statistics which it stores in a global vector of

length 2, which can be accessed by various :doc:`output commands <Howto_output>`.  The vector values calculated by this fix

are "intensive".

length 2, which can be accessed by various :doc:`output commands

<Howto_output>`.  The vector values calculated by this fix are

"intensive".

These are the 2 quantities:

#include "random_mars.h"

#include "memory.h"

#include "error.h"

#include "math_const.h"

using namespace LAMMPS_NS;

using namespace FixConst;

using namespace MathConst;

#define BIG 1.0e20

#define DELTA 16

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

      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 &&

        strcmp(style,"bond/create/angle") == 0) {

      if (iarg+3 > narg) 

          error->all(FLERR,"Illegal fix bond/create/angle 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/angle command");

      if (amin < 0 || amin > 180)

        error->all(FLERR,"Illegal fix bond/create/angle command");

      if (amax < 0 || amax > 180)

        error->all(FLERR,"Illegal fix bond/create/angle command");

      amin = (MY_PI/180.0) * amin;

      amax = (MY_PI/180.0) * amax;

      constrainflag = 1;

      iarg += 3;

    } else error->all(FLERR,"Illegal fix bond/create command");

  }

      rsq = delx*delx + dely*dely + 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;

        partner[j] = tag[i];

        distsq[j] = rsq;

      }

    }

  }

class FixBondCreate : public Fix {

 public:

  FixBondCreate(class LAMMPS *, int, char **);

  ~FixBondCreate();

  virtual ~FixBondCreate();

  int setmask();

  void init();

  void init_list(int, class NeighList *);

  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;

  void create_impropers(int);

  int dedup(int, int, tagint *);

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

  // DEBUG

  void print_bb();