forked from lijiext/lammps
added options to compute bond/local
This commit is contained in:
Steve Plimpton
parent
fe2fca4e9b
commit
25e518a4f4
|
|
@ -128,9 +128,9 @@
|
|||
<span id="index-0"></span><h1>comm_modify command</h1>
|
||||
<div class="section" id="syntax">
|
||||
<h2>Syntax</h2>
|
||||
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">comm_modify</span> <span class="n">keyword</span> <span class="n">value</span> <span class="o">...</span>
|
||||
</pre></div>
|
||||
</div>
|
||||
<pre class="literal-block">
|
||||
comm_modify keyword value ...
|
||||
</pre>
|
||||
<ul class="simple">
|
||||
<li>zero or more keyword/value pairs may be appended</li>
|
||||
<li>keyword = <em>mode</em> or <em>cutoff</em> or <em>cutoff/multi</em> or <em>group</em> or <em>vel</em></li>
|
||||
|
|
@ -147,14 +147,14 @@
|
|||
</div>
|
||||
<div class="section" id="examples">
|
||||
<h2>Examples</h2>
|
||||
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">comm_modify</span> <span class="n">mode</span> <span class="n">multi</span>
|
||||
<span class="n">comm_modify</span> <span class="n">mode</span> <span class="n">multi</span> <span class="n">group</span> <span class="n">solvent</span>
|
||||
<span class="n">comm_modift</span> <span class="n">mode</span> <span class="n">multi</span> <span class="n">cutoff</span><span class="o">/</span><span class="n">multi</span> <span class="mi">1</span> <span class="mf">10.0</span> <span class="n">cutoff</span><span class="o">/</span><span class="n">multi</span> <span class="mi">2</span><span class="o">*</span><span class="mi">4</span> <span class="mf">15.0</span>
|
||||
<span class="n">comm_modify</span> <span class="n">vel</span> <span class="n">yes</span>
|
||||
<span class="n">comm_modify</span> <span class="n">mode</span> <span class="n">single</span> <span class="n">cutoff</span> <span class="mf">5.0</span> <span class="n">vel</span> <span class="n">yes</span>
|
||||
<span class="n">comm_modify</span> <span class="n">cutoff</span><span class="o">/</span><span class="n">multi</span> <span class="o">*</span> <span class="mf">0.0</span>
|
||||
</pre></div>
|
||||
</div>
|
||||
<pre class="literal-block">
|
||||
comm_modify mode multi
|
||||
comm_modify mode multi group solvent
|
||||
comm_modift mode multi cutoff/multi 1 10.0 cutoff/multi 2*4 15.0
|
||||
comm_modify vel yes
|
||||
comm_modify mode single cutoff 5.0 vel yes
|
||||
comm_modify cutoff/multi * 0.0
|
||||
</pre>
|
||||
</div>
|
||||
<div class="section" id="description">
|
||||
<h2>Description</h2>
|
||||
|
|
|
|||
|
|
@ -135,18 +135,25 @@
|
|||
<li>ID, group-ID are documented in <a class="reference internal" href="compute.html"><span class="doc">compute</span></a> command</li>
|
||||
<li>bond/local = style name of this compute command</li>
|
||||
<li>one or more values may be appended</li>
|
||||
<li>value = <em>dist</em> or <em>eng</em> or <em>force</em></li>
|
||||
<li>value = <em>dist</em> or <em>engpot</em> or <em>force</em> or <em>engvib</em> or <em>engrot</em> or <em>engtrans</em> or <em>omega</em> or <em>velvib</em></li>
|
||||
</ul>
|
||||
<pre class="literal-block">
|
||||
<em>dist</em> = bond distance
|
||||
<em>eng</em> = bond energy
|
||||
<em>engpot</em> = bond potential energy
|
||||
<em>force</em> = bond force
|
||||
</pre>
|
||||
<pre class="literal-block">
|
||||
<em>engvib</em> = bond kinetic energy of vibration
|
||||
<em>engrot</em> = bond kinetic energy of rotation
|
||||
<em>engtrans</em> = bond kinetic energy of translation
|
||||
<em>omega</em> = magnitude of bond angular velocity
|
||||
<em>velvib</em> = vibrational velocity along the bond length
|
||||
</pre>
|
||||
</div>
|
||||
<div class="section" id="examples">
|
||||
<h2>Examples</h2>
|
||||
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">compute</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">bond</span><span class="o">/</span><span class="n">local</span> <span class="n">eng</span>
|
||||
<span class="n">compute</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">bond</span><span class="o">/</span><span class="n">local</span> <span class="n">dist</span> <span class="n">eng</span> <span class="n">force</span>
|
||||
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">compute</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">bond</span><span class="o">/</span><span class="n">local</span> <span class="n">engpot</span>
|
||||
<span class="n">compute</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">bond</span><span class="o">/</span><span class="n">local</span> <span class="n">dist</span> <span class="n">engpot</span> <span class="n">force</span>
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
|
|
@ -154,12 +161,42 @@
|
|||
<h2>Description</h2>
|
||||
<p>Define a computation that calculates properties of individual bond
|
||||
interactions. The number of datums generated, aggregated across all
|
||||
processors, equals the number of bonds in the system, modified
|
||||
by the group parameter as explained below.</p>
|
||||
<p>The value <em>dist</em> is the length of the bond.</p>
|
||||
<p>The value <em>eng</em> is the interaction energy for the bond.</p>
|
||||
<p>The value <em>force</em> is the force acting between the pair of atoms in the
|
||||
processors, equals the number of bonds in the system, modified by the
|
||||
group parameter as explained below.</p>
|
||||
<p>The value <em>dist</em> is the current length of the bond.</p>
|
||||
<p>The value <em>engpot</em> is the potential energy for the bond,
|
||||
based on the current separation of the pair of atoms in the bond.</p>
|
||||
<p>The value <em>force</em> is the magnitude of the force acting between the
|
||||
pair of atoms in the bond.</p>
|
||||
<p>The remaining properties are all computed for motion of the two atoms
|
||||
relative to the center of mass (COM) velocity of the 2 atoms in the
|
||||
bond.</p>
|
||||
<p>The value <em>engvib</em> is the vibrational kinetic energy of the two atoms
|
||||
in the bond, which is simply 1/2 m1 v1^2 + 1/2 m1 v2^2, where v1 and
|
||||
v2 are the magnitude of the velocity of the 2 atoms along the bond
|
||||
direction, after the COM velocity has been subtracted from each.</p>
|
||||
<p>The value <em>engrot</em> is the rotationsl kinetic energy of the two atoms
|
||||
in the bond, which is simply 1/2 m1 v1^2 + 1/2 m1 v2^2, where v1 and
|
||||
v2 are the magnitude of the velocity of the 2 atoms perpendicular to
|
||||
the bond direction, after the COM velocity has been subtracted from
|
||||
each.</p>
|
||||
<p>The value <em>engtrans</em> is the translational kinetic energy associated
|
||||
with the motion of the COM of the system itself, namely 1/2 (m1+m2)
|
||||
Vcm^2 where Vcm = magnitude of the velocity of the COM.</p>
|
||||
<p>Note that these 3 kinetic energy terms are simply a partitioning of
|
||||
the summed kinetic energy of the 2 atoms themselves. I.e. total KE =
|
||||
1/2 m1 v1^2 + 1/2 m2 v3^2 = engvib + engrot + engtrans, where v1,v2
|
||||
are the magnitude of the velocities of the 2 atoms, without any
|
||||
adjustment for the COM velocity.</p>
|
||||
<p>The value <em>omega</em> is the magnitude of the angular velocity of the
|
||||
two atoms around their COM position.</p>
|
||||
<p>The value <em>velvib</em> is the magnitude of the relative velocity of the
|
||||
two atoms in the bond towards each other. A negative value means the
|
||||
2 atoms are moving toward each other; a positive value means they are
|
||||
moving apart.</p>
|
||||
<p>Note that all these properties are computed for the pair of atoms in a
|
||||
bond, whether the 2 atoms represent a simple diatomic molecule, or are
|
||||
part of some larger molecule.</p>
|
||||
<p>The local data stored by this command is generated by looping over all
|
||||
the atoms owned on a processor and their bonds. A bond will only be
|
||||
included if both atoms in the bond are in the specified compute group.
|
||||
|
|
@ -179,8 +216,8 @@ command in a consistent way.</p>
|
|||
<p>Here is an example of how to do this:</p>
|
||||
<pre class="literal-block">
|
||||
compute 1 all property/local btype batom1 batom2
|
||||
compute 2 all bond/local dist eng
|
||||
dump 1 all local 1000 tmp.dump index c_1[1] c_1[2] c_1[3] c_2[1] c_2[2]
|
||||
compute 2 all bond/local dist engpot
|
||||
dump 1 all local 1000 tmp.dump index c_1[*] c_2[*]
|
||||
</pre>
|
||||
<p><strong>Output info:</strong></p>
|
||||
<p>This compute calculates a local vector or local array depending on the
|
||||
|
|
@ -191,9 +228,12 @@ local array is produced where the number of columns = the number of
|
|||
keywords. The vector or array can be accessed by any command that
|
||||
uses local values from a compute as input. See <a class="reference internal" href="Section_howto.html#howto-15"><span class="std std-ref">this section</span></a> for an overview of LAMMPS output
|
||||
options.</p>
|
||||
<p>The output for <em>dist</em> will be in distance <a class="reference internal" href="units.html"><span class="doc">units</span></a>. The
|
||||
output for <em>eng</em> will be in energy <a class="reference internal" href="units.html"><span class="doc">units</span></a>. The output for
|
||||
<em>force</em> will be in force <a class="reference internal" href="units.html"><span class="doc">units</span></a>.</p>
|
||||
<p>The output for <em>dist</em> will be in distance <a class="reference internal" href="units.html"><span class="doc">units</span></a>. The
|
||||
output for <em>velvib</em> will be in velocity <a class="reference internal" href="units.html"><span class="doc">units</span></a>. The output
|
||||
for <em>omega</em> will be in velocity/distance <a class="reference internal" href="units.html"><span class="doc">units</span></a>. The
|
||||
output for <em>engtrans</em>, <em>engvib</em>, <em>engrot</em>, and <em>engpot</em> will be in
|
||||
energy <a class="reference internal" href="units.html"><span class="doc">units</span></a>. The output for <em>force</em> will be in force
|
||||
<a class="reference internal" href="units.html"><span class="doc">units</span></a>.</p>
|
||||
</div>
|
||||
<div class="section" id="restrictions">
|
||||
<h2>Restrictions</h2>
|
||||
|
|
|
|||
|
|
@ -344,7 +344,7 @@
|
|||
</dt>
|
||||
|
||||
|
||||
<dt><a href="comm_modify.html#index-0">comm_modify</a>
|
||||
<dt><a href="comm_modify.html#index-0">comm\_modify</a>
|
||||
</dt>
|
||||
|
||||
|
||||
|
|
@ -1768,6 +1768,10 @@
|
|||
</dt>
|
||||
|
||||
|
||||
<dt><a href="pair_vashishta.html#index-0">pair\_style vashishta</a>
|
||||
</dt>
|
||||
|
||||
|
||||
<dt><a href="pair_coeff.html#index-0">pair_coeff</a>
|
||||
</dt>
|
||||
|
||||
|
|
@ -1939,12 +1943,12 @@
|
|||
<dt><a href="pair_lj_sf.html#index-0">pair_style lj/sf</a>
|
||||
</dt>
|
||||
|
||||
</dl></td>
|
||||
<td style="width: 33%" valign="top"><dl>
|
||||
|
||||
<dt><a href="pair_lj_smooth.html#index-0">pair_style lj/smooth</a>
|
||||
</dt>
|
||||
|
||||
</dl></td>
|
||||
<td style="width: 33%" valign="top"><dl>
|
||||
|
||||
<dt><a href="pair_lj_smooth_linear.html#index-0">pair_style lj/smooth/linear</a>
|
||||
</dt>
|
||||
|
|
@ -2102,10 +2106,6 @@
|
|||
</dt>
|
||||
|
||||
|
||||
<dt><a href="pair_vashishta.html#index-0">pair_style vashishta</a>
|
||||
</dt>
|
||||
|
||||
|
||||
<dt><a href="pair_yukawa.html#index-0">pair_style yukawa</a>
|
||||
</dt>
|
||||
|
||||
|
|
|
|||
|
|
@ -15,31 +15,74 @@ compute ID group-ID bond/local value1 value2 ... :pre
|
|||
ID, group-ID are documented in "compute"_compute.html command :ulb,l
|
||||
bond/local = style name of this compute command :l
|
||||
one or more values may be appended :l
|
||||
value = {dist} or {eng} or {force} :l
|
||||
value = {dist} or {engpot} or {force} or {engvib} or {engrot} or {engtrans} or {omega} or {velvib} :l
|
||||
{dist} = bond distance
|
||||
{eng} = bond energy
|
||||
{engpot} = bond potential energy
|
||||
{force} = bond force :pre
|
||||
{engvib} = bond kinetic energy of vibration
|
||||
{engrot} = bond kinetic energy of rotation
|
||||
{engtrans} = bond kinetic energy of translation
|
||||
{omega} = magnitude of bond angular velocity
|
||||
{velvib} = vibrational velocity along the bond length :pre
|
||||
:ule
|
||||
|
||||
[Examples:]
|
||||
|
||||
compute 1 all bond/local eng
|
||||
compute 1 all bond/local dist eng force :pre
|
||||
compute 1 all bond/local engpot
|
||||
compute 1 all bond/local dist engpot force :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
Define a computation that calculates properties of individual bond
|
||||
interactions. The number of datums generated, aggregated across all
|
||||
processors, equals the number of bonds in the system, modified
|
||||
by the group parameter as explained below.
|
||||
processors, equals the number of bonds in the system, modified by the
|
||||
group parameter as explained below.
|
||||
|
||||
The value {dist} is the length of the bond.
|
||||
The value {dist} is the current length of the bond.
|
||||
|
||||
The value {eng} is the interaction energy for the bond.
|
||||
The value {engpot} is the potential energy for the bond,
|
||||
based on the current separation of the pair of atoms in the bond.
|
||||
|
||||
The value {force} is the force acting between the pair of atoms in the
|
||||
The value {force} is the magnitude of the force acting between the
|
||||
pair of atoms in the bond.
|
||||
|
||||
The remaining properties are all computed for motion of the two atoms
|
||||
relative to the center of mass (COM) velocity of the 2 atoms in the
|
||||
bond.
|
||||
|
||||
The value {engvib} is the vibrational kinetic energy of the two atoms
|
||||
in the bond, which is simply 1/2 m1 v1^2 + 1/2 m1 v2^2, where v1 and
|
||||
v2 are the magnitude of the velocity of the 2 atoms along the bond
|
||||
direction, after the COM velocity has been subtracted from each.
|
||||
|
||||
The value {engrot} is the rotationsl kinetic energy of the two atoms
|
||||
in the bond, which is simply 1/2 m1 v1^2 + 1/2 m1 v2^2, where v1 and
|
||||
v2 are the magnitude of the velocity of the 2 atoms perpendicular to
|
||||
the bond direction, after the COM velocity has been subtracted from
|
||||
each.
|
||||
|
||||
The value {engtrans} is the translational kinetic energy associated
|
||||
with the motion of the COM of the system itself, namely 1/2 (m1+m2)
|
||||
Vcm^2 where Vcm = magnitude of the velocity of the COM.
|
||||
|
||||
Note that these 3 kinetic energy terms are simply a partitioning of
|
||||
the summed kinetic energy of the 2 atoms themselves. I.e. total KE =
|
||||
1/2 m1 v1^2 + 1/2 m2 v3^2 = engvib + engrot + engtrans, where v1,v2
|
||||
are the magnitude of the velocities of the 2 atoms, without any
|
||||
adjustment for the COM velocity.
|
||||
|
||||
The value {omega} is the magnitude of the angular velocity of the
|
||||
two atoms around their COM position.
|
||||
|
||||
The value {velvib} is the magnitude of the relative velocity of the
|
||||
two atoms in the bond towards each other. A negative value means the
|
||||
2 atoms are moving toward each other; a positive value means they are
|
||||
moving apart.
|
||||
|
||||
Note that all these properties are computed for the pair of atoms in a
|
||||
bond, whether the 2 atoms represent a simple diatomic molecule, or are
|
||||
part of some larger molecule.
|
||||
|
||||
The local data stored by this command is generated by looping over all
|
||||
the atoms owned on a processor and their bonds. A bond will only be
|
||||
included if both atoms in the bond are in the specified compute group.
|
||||
|
|
@ -62,8 +105,8 @@ command in a consistent way.
|
|||
Here is an example of how to do this:
|
||||
|
||||
compute 1 all property/local btype batom1 batom2
|
||||
compute 2 all bond/local dist eng
|
||||
dump 1 all local 1000 tmp.dump index c_1\[1\] c_1\[2\] c_1\[3\] c_2\[1\] c_2\[2\] :pre
|
||||
compute 2 all bond/local dist engpot
|
||||
dump 1 all local 1000 tmp.dump index c_1\[*\] c_2\[*\] :pre
|
||||
|
||||
[Output info:]
|
||||
|
||||
|
|
@ -77,9 +120,12 @@ uses local values from a compute as input. See "this
|
|||
section"_Section_howto.html#howto_15 for an overview of LAMMPS output
|
||||
options.
|
||||
|
||||
The output for {dist} will be in distance "units"_units.html. The
|
||||
output for {eng} will be in energy "units"_units.html. The output for
|
||||
{force} will be in force "units"_units.html.
|
||||
The output for {dist} will be in distance "units"_units.html. The
|
||||
output for {velvib} will be in velocity "units"_units.html. The output
|
||||
for {omega} will be in velocity/distance "units"_units.html. The
|
||||
output for {engtrans}, {engvib}, {engrot}, and {engpot} will be in
|
||||
energy "units"_units.html. The output for {force} will be in force
|
||||
"units"_units.html.
|
||||
|
||||
[Restrictions:] none
|
||||
|
||||
|
|
|
|||
|
|
@ -21,14 +21,17 @@
|
|||
#include "domain.h"
|
||||
#include "force.h"
|
||||
#include "bond.h"
|
||||
#include "math_extra.h"
|
||||
#include "comm.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define DELTA 10000
|
||||
#define EPSILON 1.0e-12
|
||||
|
||||
enum{DIST,ENG,FORCE};
|
||||
enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
|
|
@ -42,6 +45,8 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
|
|||
error->all(FLERR,"Compute bond/local used when bonds are not allowed");
|
||||
|
||||
local_flag = 1;
|
||||
comm_forward = 3;
|
||||
|
||||
nvalues = narg - 3;
|
||||
if (nvalues == 1) size_local_cols = 0;
|
||||
else size_local_cols = nvalues;
|
||||
|
|
@ -51,16 +56,26 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
|
|||
nvalues = 0;
|
||||
for (int iarg = 3; iarg < narg; iarg++) {
|
||||
if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST;
|
||||
else if (strcmp(arg[iarg],"eng") == 0) bstyle[nvalues++] = ENG;
|
||||
else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT;
|
||||
else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE;
|
||||
else if (strcmp(arg[iarg],"engvib") == 0) bstyle[nvalues++] = ENGVIB;
|
||||
else if (strcmp(arg[iarg],"engrot") == 0) bstyle[nvalues++] = ENGROT;
|
||||
else if (strcmp(arg[iarg],"engtrans") == 0) bstyle[nvalues++] = ENGTRANS;
|
||||
else if (strcmp(arg[iarg],"omega") == 0) bstyle[nvalues++] = OMEGA;
|
||||
else if (strcmp(arg[iarg],"velvib") == 0) bstyle[nvalues++] = VELVIB;
|
||||
else error->all(FLERR,"Invalid keyword in compute bond/local command");
|
||||
}
|
||||
|
||||
// set singleflag if need to call bond->single()
|
||||
// set velflag if compute any quantities based on velocities
|
||||
|
||||
singleflag = 0;
|
||||
for (int i = 0; i < nvalues; i++)
|
||||
if (bstyle[i] != DIST) singleflag = 1;
|
||||
ghostvelflag = 0;
|
||||
for (int i = 0; i < nvalues; i++) {
|
||||
if (bstyle[i] == ENGPOT || bstyle[i] == FORCE) singleflag = 1;
|
||||
if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS ||
|
||||
bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1;
|
||||
}
|
||||
|
||||
nmax = 0;
|
||||
}
|
||||
|
|
@ -81,9 +96,17 @@ void ComputeBondLocal::init()
|
|||
if (force->bond == NULL)
|
||||
error->all(FLERR,"No bond style is defined for compute bond/local");
|
||||
|
||||
// set ghostvelflag if need to acquire ghost atom velocities
|
||||
|
||||
if (velflag && !comm->ghost_velocity) ghostvelflag = 1;
|
||||
else ghostvelflag = 0;
|
||||
|
||||
// do initial memory allocation so that memory_usage() is correct
|
||||
|
||||
initflag = 1;
|
||||
ncount = compute_bonds(0);
|
||||
initflag = 0;
|
||||
|
||||
if (ncount > nmax) reallocate(ncount);
|
||||
size_local_rows = ncount;
|
||||
}
|
||||
|
|
@ -117,10 +140,22 @@ int ComputeBondLocal::compute_bonds(int flag)
|
|||
{
|
||||
int i,m,n,nb,atom1,atom2,imol,iatom,btype;
|
||||
tagint tagprev;
|
||||
double delx,dely,delz,rsq;
|
||||
double dx,dy,dz,rsq;
|
||||
double mass1,mass2,masstotal,invmasstotal;
|
||||
double xcm[3],vcm[3];
|
||||
double delr1[3],delr2[3],delv1[3],delv2[3];
|
||||
double r12[3],vpar1,vpar2;
|
||||
double vvib,vrotsq;
|
||||
double inertia,omegasq;
|
||||
double mvv2e;
|
||||
double engpot,engtrans,engvib,engrot,engtot,fbond;
|
||||
double *ptr;
|
||||
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
int *type = atom->type;
|
||||
double *mass = atom->mass;
|
||||
double *rmass = atom->rmass;
|
||||
tagint *tag = atom->tag;
|
||||
int *num_bond = atom->num_bond;
|
||||
tagint **bond_atom = atom->bond_atom;
|
||||
|
|
@ -136,7 +171,12 @@ int ComputeBondLocal::compute_bonds(int flag)
|
|||
int molecular = atom->molecular;
|
||||
|
||||
Bond *bond = force->bond;
|
||||
double eng,fbond;
|
||||
|
||||
// communicate ghost velocities if needed
|
||||
|
||||
if (ghostvelflag && !initflag) comm->forward_comm_compute(this);
|
||||
|
||||
// loop over all atoms and their bonds
|
||||
|
||||
m = n = 0;
|
||||
for (atom1 = 0; atom1 < nlocal; atom1++) {
|
||||
|
|
@ -164,16 +204,93 @@ int ComputeBondLocal::compute_bonds(int flag)
|
|||
if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
|
||||
if (btype == 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 = delx*delx + dely*dely + delz*delz;
|
||||
if (!flag) {
|
||||
m++;
|
||||
continue;
|
||||
}
|
||||
|
||||
if (singleflag && (btype > 0))
|
||||
eng = bond->single(btype,rsq,atom1,atom2,fbond);
|
||||
else eng = fbond = 0.0;
|
||||
dx = x[atom1][0] - x[atom2][0];
|
||||
dy = x[atom1][1] - x[atom2][1];
|
||||
dz = x[atom1][2] - x[atom2][2];
|
||||
domain->minimum_image(dx,dy,dz);
|
||||
rsq = dx*dx + dy*dy + dz*dz;
|
||||
|
||||
if (btype == 0) {
|
||||
engpot = fbond = 0.0;
|
||||
engvib = engrot = engtrans = omegasq = vvib = 0.0;
|
||||
} else {
|
||||
|
||||
if (singleflag) engpot = bond->single(btype,rsq,atom1,atom2,fbond);
|
||||
|
||||
if (velflag) {
|
||||
if (rmass) {
|
||||
mass1 = rmass[atom1];
|
||||
mass2 = rmass[atom2];
|
||||
}
|
||||
else {
|
||||
mass1 = mass[type[atom1]];
|
||||
mass2 = mass[type[atom2]];
|
||||
}
|
||||
masstotal = mass1+mass2;
|
||||
invmasstotal = 1.0 / (masstotal);
|
||||
xcm[0] = (mass1*x[atom1][0] + mass2*x[atom2][0]) * invmasstotal;
|
||||
xcm[1] = (mass1*x[atom1][1] + mass2*x[atom2][1]) * invmasstotal;
|
||||
xcm[2] = (mass1*x[atom1][2] + mass2*x[atom2][2]) * invmasstotal;
|
||||
vcm[0] = (mass1*v[atom1][0] + mass2*v[atom2][0]) * invmasstotal;
|
||||
vcm[1] = (mass1*v[atom1][1] + mass2*v[atom2][1]) * invmasstotal;
|
||||
vcm[2] = (mass1*v[atom1][2] + mass2*v[atom2][2]) * invmasstotal;
|
||||
|
||||
engtrans = 0.5 * masstotal * MathExtra::lensq3(vcm);
|
||||
|
||||
// r12 = unit bond vector from atom1 to atom2
|
||||
|
||||
MathExtra::sub3(x[atom2],x[atom1],r12);
|
||||
MathExtra::norm3(r12);
|
||||
|
||||
// delr = vector from COM to each atom
|
||||
// delv = velocity of each atom relative to COM
|
||||
|
||||
MathExtra::sub3(x[atom1],xcm,delr1);
|
||||
MathExtra::sub3(x[atom2],xcm,delr2);
|
||||
MathExtra::sub3(v[atom1],vcm,delv1);
|
||||
MathExtra::sub3(v[atom2],vcm,delv2);
|
||||
|
||||
// vpar = component of delv parallel to bond vector
|
||||
|
||||
vpar1 = MathExtra::dot3(delv1,r12);
|
||||
vpar2 = MathExtra::dot3(delv2,r12);
|
||||
engvib = 0.5 * (mass1*vpar1*vpar1 + mass2*vpar2*vpar2);
|
||||
|
||||
// vvib = relative velocity of 2 atoms along bond direction
|
||||
// vvib < 0 for 2 atoms moving towards each other
|
||||
// vvib > 0 for 2 atoms moving apart
|
||||
|
||||
vvib = vpar2 - vpar1;
|
||||
|
||||
// vrotsq = tangential speed squared of atom1 only
|
||||
// omegasq = omega squared, and is the same for atom1 and atom2
|
||||
|
||||
inertia = mass1*MathExtra::lensq3(delr1) +
|
||||
mass2*MathExtra::lensq3(delr2);
|
||||
vrotsq = MathExtra::lensq3(delv1) - vpar1*vpar1;
|
||||
omegasq = vrotsq / MathExtra::lensq3(delr1);
|
||||
|
||||
engrot = 0.5 * inertia * omegasq;
|
||||
|
||||
// sanity check: engtotal = engtrans + engvib + engrot
|
||||
|
||||
//engtot = 0.5 * (mass1*MathExtra::lensq3(v[atom1]) +
|
||||
// mass2*MathExtra::lensq3(v[atom2]));
|
||||
//if (fabs(engtot-engtrans-engvib-engrot) > EPSILON)
|
||||
// error->one(FLERR,"Sanity check on 3 energy components failed");
|
||||
|
||||
// scale energies by units
|
||||
|
||||
mvv2e = force->mvv2e;
|
||||
engtrans *= mvv2e;
|
||||
engvib *= mvv2e;
|
||||
engrot *= mvv2e;
|
||||
}
|
||||
|
||||
if (nvalues == 1) ptr = &vector[m];
|
||||
else ptr = array[m];
|
||||
|
|
@ -183,12 +300,27 @@ int ComputeBondLocal::compute_bonds(int flag)
|
|||
case DIST:
|
||||
ptr[n] = sqrt(rsq);
|
||||
break;
|
||||
case ENG:
|
||||
ptr[n] = eng;
|
||||
case ENGPOT:
|
||||
ptr[n] = engpot;
|
||||
break;
|
||||
case FORCE:
|
||||
ptr[n] = sqrt(rsq)*fbond;
|
||||
break;
|
||||
case ENGVIB:
|
||||
ptr[n] = engvib;
|
||||
break;
|
||||
case ENGROT:
|
||||
ptr[n] = engrot;
|
||||
break;
|
||||
case ENGTRANS:
|
||||
ptr[n] = engtrans;
|
||||
break;
|
||||
case OMEGA:
|
||||
ptr[n] = sqrt(omegasq);
|
||||
break;
|
||||
case VELVIB:
|
||||
ptr[n] = vvib;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -202,6 +334,43 @@ int ComputeBondLocal::compute_bonds(int flag)
|
|||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf,
|
||||
int pbc_flag, int *pbc)
|
||||
{
|
||||
int i,j,m;
|
||||
|
||||
double **v = atom->v;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
buf[m++] = v[j][0];
|
||||
buf[m++] = v[j][1];
|
||||
buf[m++] = v[j][2];
|
||||
}
|
||||
|
||||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeBondLocal::unpack_forward_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i,m,last;
|
||||
|
||||
double **v = atom->v;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++) {
|
||||
v[i][0] = buf[m++];
|
||||
v[i][1] = buf[m++];
|
||||
v[i][2] = buf[m++];
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeBondLocal::reallocate(int n)
|
||||
{
|
||||
// grow vector or array and indices array
|
||||
|
|
|
|||
|
|
@ -30,13 +30,15 @@ class ComputeBondLocal : public Compute {
|
|||
~ComputeBondLocal();
|
||||
void init();
|
||||
void compute_local();
|
||||
int pack_forward_comm(int, int *, double *, int, int *);
|
||||
void unpack_forward_comm(int, int, double *);
|
||||
double memory_usage();
|
||||
|
||||
private:
|
||||
int nvalues;
|
||||
int ncount;
|
||||
int *bstyle;
|
||||
int singleflag;
|
||||
int singleflag,velflag,ghostvelflag,initflag;
|
||||
|
||||
int nmax;
|
||||
|
||||
|
|
|
|||
|
|
@ -21,6 +21,7 @@
|
|||
#include "compute.h"
|
||||
#include "domain.h"
|
||||
#include "update.h"
|
||||
#include "input.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
|
|
@ -45,7 +46,18 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) :
|
|||
nevery = force->inumeric(FLERR,arg[3]);
|
||||
if (nevery <= 0) error->all(FLERR,"Illegal dump local command");
|
||||
|
||||
size_one = nfield = narg-5;
|
||||
nfield = narg - 5;
|
||||
|
||||
// expand args if any have wildcard character "*"
|
||||
|
||||
int expand = 0;
|
||||
char **earg;
|
||||
nfield = input->expand_args(nfield,&arg[5],1,earg);
|
||||
|
||||
if (earg != &arg[5]) expand = 1;
|
||||
|
||||
// allocate field vectors
|
||||
|
||||
pack_choice = new FnPtrPack[nfield];
|
||||
vtype = new int[nfield];
|
||||
|
||||
|
|
@ -67,7 +79,8 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) :
|
|||
|
||||
// process attributes
|
||||
|
||||
parse_fields(narg,arg);
|
||||
parse_fields(nfield,earg);
|
||||
size_one = nfield;
|
||||
|
||||
// setup format strings
|
||||
|
||||
|
|
@ -88,11 +101,11 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) :
|
|||
// setup column string
|
||||
|
||||
int n = 0;
|
||||
for (int iarg = 5; iarg < narg; iarg++) n += strlen(arg[iarg]) + 2;
|
||||
for (int iarg = 0; iarg < nfield; iarg++) n += strlen(earg[iarg]) + 2;
|
||||
columns = new char[n];
|
||||
columns[0] = '\0';
|
||||
for (int iarg = 5; iarg < narg; iarg++) {
|
||||
strcat(columns,arg[iarg]);
|
||||
for (int iarg = 0; iarg < nfield; iarg++) {
|
||||
strcat(columns,earg[iarg]);
|
||||
strcat(columns," ");
|
||||
}
|
||||
|
||||
|
|
@ -102,6 +115,13 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) :
|
|||
n = strlen(str) + 1;
|
||||
label = new char[n];
|
||||
strcpy(label,str);
|
||||
|
||||
// if wildcard expansion occurred, free earg memory from exapnd_args()
|
||||
|
||||
if (expand) {
|
||||
for (int i = 0; i < nfield; i++) delete [] earg[i];
|
||||
memory->sfree(earg);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -376,8 +396,8 @@ void DumpLocal::parse_fields(int narg, char **arg)
|
|||
// customize by adding to if statement
|
||||
|
||||
int i;
|
||||
for (int iarg = 5; iarg < narg; iarg++) {
|
||||
i = iarg-5;
|
||||
for (int iarg = 0; iarg < narg; iarg++) {
|
||||
i = iarg;
|
||||
|
||||
if (strcmp(arg[iarg],"index") == 0) {
|
||||
pack_choice[i] = &DumpLocal::pack_index;
|
||||
|
|
|
|||
|
|
@ -36,7 +36,8 @@ enum{X,V,F,COMPUTE,FIX,VARIABLE};
|
|||
|
||||
FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg),
|
||||
nvalues(0), which(NULL), argindex(NULL), value2index(NULL), ids(NULL), array(NULL)
|
||||
nvalues(0), which(NULL), argindex(NULL), value2index(NULL),
|
||||
ids(NULL), array(NULL)
|
||||
{
|
||||
if (narg < 7) error->all(FLERR,"Illegal fix ave/atom command");
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue