Homec4science
Diffusion lammps (lammm-master-deprecated) Browse

lammps/tools/moltemplate/src/nbody_Angles.py
d510bb54ef32lammm-master-deprecated

Last Change
Raw File

nbody_Angles.py
View Options

from nbody_graph_search import Ugraph
# This file defines how 3-body angle interactions are generated by moltemplate
# by default. It can be overridden by supplying your own custom file.
# To find 3-body "angle" interactions, we would use this subgraph:
#
#
# *---*---* => 1st bond connects atoms 0 and 1
# 0 1 2 2nd bond connects atoms 1 and 2
#
bond_pattern = Ugraph([(0,1), (1,2)])
# (Ugraph atom indices begin at 0, not 1)
# The next function eliminates the redundancy between 0-1-2 and 2-1-0:
def canonical_order(match):
"""
Before defining a new interaction, we must check to see if an
interaction between these same 3 atoms has already been created
(perhaps listed in a different, but equivalent order).
If we don't check for this this, we will create many unnecessary redundant
interactions (which can slow down he simulation).
To avoid this, I define a "canonical_order" function which sorts the atoms
and bonds in a way which is consistent with the symmetry of the interaction
being generated... Later the re-ordered list of atom and bond ids will be
tested against the list of atom/bond ids in the matches-found-so-far,
before it is added to the list of interactions found so far. Note that
the energy of an angle interaction is a function of the angle between.
three consecutively bonded atoms (referred to here as: 0,1,2).
This angle does not change when swapping the atoms at either end (0 and 2).
So it does not make sense to define a separate 3-body angle
interaction between atoms 0,1,2 AS WELL AS an interaction between 2,1,0.
So we sort the atoms and bonds so that the first atom has a always has
a lower atomID than the third atom. (Later we will check to see if we
have already defined an interaction between these 3 atoms. If not then
we create a new one.)
"""
# match[0][0:2] contains the ID numbers for the 3 atoms in the match
atom0 = match[0][0]
atom1 = match[0][1]
atom2 = match[0][2]
# match[1][0:1] contains the ID numbers for the 2 bonds
bond0 = match[1][0]
bond1 = match[1][1]
if atom0 < atom2:
#return ((atom0, atom1, atom2), (bond0, bond1)) same thing as:
return match
else:
return ((atom2, atom1, atom0), (bond1, bond0))
c4science · Help