Commit 40088f55 authored by Oliver Henrich's avatar Oliver Henrich
Browse files

Revert "Added fix bond/create/angle code and docu"

This reverts commit 16eab647.
parent f0a32561
Loading
Loading
Loading
Loading
+0 −1
Original line number Diff line number Diff line
@@ -42,7 +42,6 @@ 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>`
+5 −23
Original line number Diff line number Diff line
@@ -2,8 +2,6 @@

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

Syntax
""""""
@@ -19,7 +17,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* or *aconstrain*
* keyword = *iparam* or *jparam* or *prob* or *atype* or *dtype* or *itype*

  .. parsed-literal::

@@ -38,9 +36,6 @@ 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
""""""""
@@ -50,7 +45,6 @@ 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
"""""""""""
@@ -116,16 +110,7 @@ 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.

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*.
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
@@ -252,11 +237,8 @@ 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.

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

built with that package.  See the :doc:`Build package <Build_package>`
doc page for more info.

Related commands
""""""""""""""""
+0 −32
Original line number Diff line number Diff line
@@ -79,11 +79,6 @@ 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) {
@@ -125,20 +120,6 @@ 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");
  }

@@ -453,11 +434,6 @@ void FixBondCreate::post_integrate()
      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;
@@ -466,7 +442,6 @@ void FixBondCreate::post_integrate()
        partner[j] = tag[i];
        distsq[j] = rsq;
      }
  
    }
  }

@@ -1431,13 +1406,6 @@ double FixBondCreate::memory_usage()

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

int FixBondCreate::constrain(int i, int j, double amin, double amax)
{
  return 1;
}

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

void FixBondCreate::print_bb()
{
  for (int i = 0; i < atom->nlocal; i++) {
+0 −5
Original line number Diff line number Diff line
@@ -52,12 +52,9 @@ class FixBondCreate : public Fix {
  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;

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

  virtual int constrain(int, int, double, double);

  // DEBUG

  void print_bb();

src/MC/fix_bond_create_angle.cpp

deleted100644 → 0
+0 −151
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)
{

}

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

FixBondCreateAngle::~FixBondCreateAngle()
{

}

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

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
    v1x = x[i][0] - x[j][0];
    v1y = x[i][1] - x[j][1];
    v1z = x[i][2] - x[j][2];

    // calculate second vector between i and its first neighbor
    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((v1x*v2x + v1y*v2y + v1z*v2z)/
        (sqrt(v1x*v1x + v1y*v1y + v1z*v1z)*sqrt(v2x*v2x + v2y*v2y + v2z*v2z)));
    // set flag if the angle constraint is met
    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
    v1x = x[j][0] - x[i][0];
    v1y = x[j][1] - x[i][1];
    v1z = x[j][2] - x[i][2];
    // calculate second vector between j and its first neighbor
    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
    angle1 = acos((v1x*v2x + v1y*v2y + v1z*v2z)/
        (sqrt(v1x*v1x + v1y*v1y + v1z*v1z)*sqrt(v2x*v2x + v2y*v2y + v2z*v2z)));
    // set flag if angle constraint is met
    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
    v1x = x[i][0] - x[j][0];
    v1y = x[i][1] - x[j][1];
    v1z = x[i][2] - x[j][2];

    // calculate second vector between i and its first neighbor
    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((v1x*v2x + v1y*v2y + v1z*v2z)/
        (sqrt(v1x*v1x + v1y*v1y + v1z*v1z)*sqrt(v2x*v2x + v2y*v2y + v2z*v2z)));

    // Calculate 2nd angle
    // calculate first vector between i and j
    v1x = x[j][0] - x[i][0];
    v1y = x[j][1] - x[i][1];
    v1z = x[j][2] - x[i][2];

    // calculate second vector between i and its first neighbor
    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((v1x*v2x + v1y*v2y + v1z*v2z)/
        (sqrt(v1x*v1x + v1y*v1y + v1z*v1z)*sqrt(v2x*v2x + v2y*v2y + v2z*v2z)));
    // set flag if the angle constraint is met

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

  }

  return flag;

}
Loading