jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' into restart-consistency

This commit is contained in:
Axel Kohlmeyer 2020-03-19 00:54:04 -04:00
parent 3e36acc533 17dd794514
commit b40bc3993d
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
33 changed files with 86526 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_fix.rst
View File

@ -119,6 +119,7 @@ OPT.
* :doc:`npt/eff <fix_nh_eff>`
* :doc:`npt/sphere (o) <fix_npt_sphere>`
* :doc:`npt/uef <fix_nh_uef>`
* :doc:`numdiff <fix_numdiff>`
* :doc:`nve (iko) <fix_nve>`
* :doc:`nve/asphere (i) <fix_nve_asphere>`
* :doc:`nve/asphere/noforce <fix_nve_asphere_noforce>`

2
doc/src/Commands_pair.rst
View File

@ -77,6 +77,8 @@ OPT.
* :doc:`coul/long/cs (g) <pair_cs>`
* :doc:`coul/long/soft (o) <pair_fep_soft>`
* :doc:`coul/msm (o) <pair_coul>`
* :doc:`coul/slater/cut <pair_coul_slater>`
* :doc:`coul/slater/long <pair_coul_slater>`
* :doc:`coul/shield <pair_coul_shield>`
* :doc:`coul/streitz <pair_coul>`
* :doc:`coul/wolf (ko) <pair_coul>`

20
doc/src/Errors_messages.rst
View File

@ -1936,6 +1936,9 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
*Compute ID for fix ave/time does not exist*
Self-explanatory.
*Compute ID for fix numdiff does not exist*
Self-explanatory.
*Compute ID for fix store/state does not exist*
Self-explanatory.
@ -3783,6 +3786,14 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
Self-explanatory. The change in the box tilt is too extreme
on a short timescale.
*Fix numdiff requires an atom map, see atom_modify*
Self-explanatory. Efficient loop over all atoms for numerical
difference requires an atom map.
*Fix numdiff requires consecutive atom IDs*
Self-explanatory. Efficient loop over all atoms for numerical
difference requires consecutive atom IDs.
*Fix nve/asphere requires extended particles*
This fix can only be used for particles with a shape setting.
@ -4733,6 +4744,12 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
*Invalid Masses section in molecule file*
Self-explanatory.
*Invalid molecule ID in molecule file*
Molecule ID must be a non-zero positive integer.
*Invalid Molecules section in molecule file*
Self-explanatory.
*Invalid REAX atom type*
There is a mis-match between LAMMPS atom types and the elements
listed in the ReaxFF force field file.
@ -4799,6 +4816,9 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
Atom IDs must be positive integers and within range of defined
atoms.
*Invalid atom ID in Fragments section of molecule file*
Self-explanatory.
*Invalid atom ID in Impropers section of data file*
Atom IDs must be positive integers and within range of defined
atoms.

2
doc/src/dynamical_matrix.rst
View File

@ -70,7 +70,7 @@ See the :doc:`Build package <Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`fix phonon <fix_phonon>`
:doc:`fix phonon <fix_phonon>`, :doc:`fix numdiff <fix_numdiff>`,
:doc:`compute hma <compute_hma>` uses an analytic formulation of the
Hessian provided by a pair_style's Pair::single_hessian() function,

1
doc/src/fix.rst
View File

@ -262,6 +262,7 @@ accelerated styles exist.
* :doc:`npt/eff <fix_nh_eff>` - NPT for nuclei and electrons in the electron force field model
* :doc:`npt/sphere <fix_npt_sphere>` - NPT for spherical particles
* :doc:`npt/uef <fix_nh_uef>` - NPT style time integration with diagonal flow
* :doc:`numdiff <fix_numdiff>` - compute derivatives of per-atom data from finite differences
* :doc:`nve <fix_nve>` - constant NVE time integration
* :doc:`nve/asphere <fix_nve_asphere>` - NVE for aspherical particles
* :doc:`nve/asphere/noforce <fix_nve_asphere_noforce>` - NVE for aspherical particles without forces

110
doc/src/fix_numdiff.rst Normal file
View File

@ -0,0 +1,110 @@
.. index:: fix numdiff
fix numdiff command
====================
Syntax
""""""
.. parsed-literal::
fix ID group-ID numdiff Nevery delta
* ID, group-ID are documented in :doc:`fix <fix>` command
* numdiff = style name of this fix command
* Nevery = calculate force by finite difference every this many timesteps
* delta = finite difference displacement length (distance units)
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all numdiff 1 0.0001
fix 1 all numdiff 10 1e-6
fix 1 all numdiff 100 0.01
Description
"""""""""""
Calculate forces through finite difference calculations of energy
versus position. These forces can be compared to analytic forces
computed by pair styles, bond styles, etc. E.g. for debugging
purposes.
The group specified with the command means only atoms within the group
have their averages computed. Results are set to 0.0 for atoms not in
the group.
This fix performs a loop over all atoms (in the group). For each atom
and each component of force it adds *delta* to the position, and
computes the new energy of the entire system. It then subtracts
*delta* from the original position and again computes the new energy
of the system. It then restores the original position. That
component of force is calculated as the difference in energy divided
by two times *delta*.
.. note::
It is important to choose a suitable value for delta, the magnitude of
atom displacements that are used to generate finite difference
approximations to the exact forces. For typical systems, a value in
the range of 1 part in 1e4 to 1e5 of the typical separation distance
between atoms in the liquid or solid state will be sufficient.
However, the best value will depend on a multitude of factors
including the stiffness of the interatomic potential, the thermodynamic
state of the material being probed, and so on. The only way to be sure
that you have made a good choice is to do a sensitivity study on a
representative atomic configuration, sweeping over a wide range of
values of delta. If delta is too small, the output forces will vary
erratically due to truncation effects. If delta is increased beyond a
certain point, the output forces will start to vary smoothly with
delta, due to growing contributions from higher order derivatives. In
between these two limits, the numerical force values should be largely
independent of delta.
.. note::
The cost of each energy evaluation is essentially the cost of an MD
timestep. This invoking this fix once has a cost of 2N timesteps,
where N is the total number of atoms in the system (assuming all atoms
are included in the group). So this fix can be very expensive to use
for large systems.
----------
The *Nevery* argument specifies on what timesteps the force will
be used calculated by finite difference.
The *delta* argument specifies the positional displacement each
atom will undergo.
----------
**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.
This fix produces a per-atom array which can be accessed by various
:doc:`output commands <Howto_output>`, which stores the components of
the force on each atom as calculated by finite difference. The
per-atom values can only be accessed on timesteps that are multiples
of *Nevery* since that is when the finite difference forces are
calculated.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is invoked during :doc:`energy
minimization <minimize>`.
Restrictions
""""""""""""
none
Related commands
""""""""""""""""
:doc:`dynamical_matrix <dynamical_matrix>`,
**Default:** none

24
doc/src/molecule.rst
View File

@ -60,6 +60,7 @@ templates include:
* :doc:`fix rigid/small <fix_rigid>`
* :doc:`fix shake <fix_shake>`
* :doc:`fix gcmc <fix_gcmc>`
* :doc:`fix bond/react <fix_bond_react>`
* :doc:`create_atoms <create_atoms>`
* :doc:`atom_style template <atom_style>`
@ -144,6 +145,7 @@ appear if the value(s) are different than the default.
* Na *angles* = # of angles Na in molecule, default = 0
* Nd *dihedrals* = # of dihedrals Nd in molecule, default = 0
* Ni *impropers* = # of impropers Ni in molecule, default = 0
* Nf *fragments* = # of fragments in molecule, default = 0
* Mtotal *mass* = total mass of molecule
* Xc Yc Zc *com* = coordinates of center-of-mass of molecule
* Ixx Iyy Izz Ixy Ixz Iyz *inertia* = 6 components of inertia tensor of molecule
@ -166,7 +168,7 @@ internally.
These are the allowed section keywords for the body of the file.
* *Coords, Types, Charges, Diameters, Masses* = atom-property sections
* *Coords, Types, Molecules, Fragments, Charges, Diameters, Masses* = atom-property sections
* *Bonds, Angles, Dihedrals, Impropers* = molecular topology sections
* *Special Bond Counts, Special Bonds* = special neighbor info
* *Shake Flags, Shake Atoms, Shake Bond Types* = SHAKE info
@ -223,6 +225,26 @@ listed in order from 1 to Nlines, but LAMMPS does not check for this.
----------
*Molecules* section:
* one line per atom
* line syntax: ID molecule-ID
* molecule-ID = molecule ID of atom
----------
*Fragments* section:
* one line per fragment
* line syntax: ID a b c d ...
* a,b,c,d,... = IDs of atoms in fragment
The ID of a fragment can only contain alphanumeric characters and
underscores. The atom IDs should be values from 1 to Natoms, where
Natoms = # of atoms in the molecule.
----------
*Charges* section:
* one line per atom

121
doc/src/pair_coul_slater.rst Normal file
View File

@ -0,0 +1,121 @@
.. index:: pair_style coul/slater
pair_style coul/slater/cut command
==================================
pair_style coul/slater/long command
===================================
Syntax
""""""
.. code-block:: LAMMPS
pair_style coul/slater/cut lamda cutoff
pair_style coul/slater/long lamda cutoff
lamda = decay length of the charge (distance units)
cutoff = cutoff (distance units)
Examples
""""""""
.. code-block:: LAMMPS
pair_style coul/slater/cut 1.0 3.5
pair_coeff * *
pair_coeff 2 2 2.5
pair_style coul/slater/long 1.0 12.0
pair_coeff * *
pair_coeff 1 1 5.0
Description
"""""""""""
Styles *coul/slater* compute electrostatic interactions in mesoscopic models
which employ potentials without explicit excluded-volume interactions.
The goal is to prevent artificial ionic pair formation by including a charge
distribution in the Coulomb potential, following the formulation of
:ref:`(Melchor) <Melchor>`:
.. math::
E = \frac{Cq_iq_j}{\epsilon r} \left( 1- \left( 1 + \frac{r_{ij}}{\lambda} exp\left( -2r_{ij}/\lambda \right) \right) \right) \qquad r < r_c
where :math:`r_c` is the cutoff distance and :math:`\lambda` is the decay length of the charge.
C is the same Coulomb conversion factor as in the pair\_styles coul/cut and coul/long. In this way the Coulomb
interaction between ions is corrected at small distances r.
For the *coul/slater/cut* style, the potential energy for distances larger than the cutoff is zero,
while for the *coul/slater/long*, the long-range interactions are computed either by the Ewald or the PPPM technique.
Phenomena that can be captured at a mesoscopic level using this type of electrostatic
interactions include the formation of polyelectrolyte-surfactant aggregates,
charge stabilization of colloidal suspensions, and the formation of
complexes driven by charged species in biological systems. :ref:`(Vaiwala) <Vaiwala>`.
The cutoff distance is optional. If it is not used,
the default global value specified in the pair_style command is used.
For each pair of atom types, a specific cutoff distance can be defined via the :doc:`pair_coeff <pair_coeff>` command as in the example
above, or in the data file or restart files read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands:
* :math:`r_c` (distance units)
The global decay length of the charge (:math:`\lambda`) specified in the pair\_style command is used for all pairs.
----------
**Mixing, shift, table, tail correction, restart, rRESPA info**\ :
For atom type pairs I,J and I != J, the cutoff distance for the
*coul/slater* styles can be mixed. The default mix value is *geometric*\ .
See the "pair\_modify" command for details.
The :doc:`pair_modify <pair_modify>` shift and table options are not relevant
for these pair styles.
These pair styles do not support the :doc:`pair_modify <pair_modify>`
tail option for adding long-range tail corrections to energy and
pressure.
These pair styles write their information to :doc:`binary restart files <restart>`, so pair\_style and pair\_coeff commands do not need
to be specified in an input script that reads a restart file.
This pair style can only be used via the *pair* keyword of the
:doc:`run_style respa <run_style>` command. It does not support the
*inner*\ , *middle*\ , *outer* keywords.
Restrictions
""""""""""""
The *coul/slater/long* style requires the long-range solvers included in the KSPACE package.
These styles are part of the "USER-MISC" package. They are only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`, :doc:`pair_style, hybrid/overlay <pair_hybrid>`, :doc:`kspace_style <kspace_style>`
**Default:** none
----------
.. _Melchor:
**(Melchor)** Gonzalez-Melchor, Mayoral, Velázquez, and Alejandre, J Chem Phys, 125, 224107 (2006).
.. _Vaiwala:
**(Vaiwala)** Vaiwala, Jadhav, and Thaokar, J Chem Phys, 146, 124904 (2017).

10
doc/utils/sphinx-config/false_positives.txt
View File

@ -1027,6 +1027,7 @@ Gmask
gneb
GNEB
Goldfarb
Gonzalez-Melchor
googlemail
Gordan
GPa
@ -1290,6 +1291,7 @@ Izvekov
izz
Izz
Jacobsen
Jadhav
jagreat
Jalalvand
james
@ -1663,6 +1665,7 @@ maxwell
Maxwellian
maxX
Mayergoyz
Mayoral
mbt
Mbytes
MBytes
@ -1691,6 +1694,7 @@ mediumvioletred
Mees
Mehl
Mei
Melchor
Meloni
Melrose
Mem
@ -2049,6 +2053,7 @@ nucleotides
num
numa
numactl
numdiff
numericalfreedom
numerics
numpy
@ -2260,6 +2265,7 @@ polyA
polybond
polydisperse
polydispersity
polyelectrolyte
polyhedra
popen
Popov
@ -2769,6 +2775,7 @@ superset
supersphere
Supinski
Surblys
surfactant
surfactants
Suter
Sutmann
@ -2831,6 +2838,7 @@ tfmc
tfMC
th
Thakkar
Thaokar
thb
thei
Theodorou
@ -3025,6 +3033,7 @@ uwo
Uzdin
vacf
Vaid
Vaiwala
valent
Valeriu
valgrind
@ -3053,6 +3062,7 @@ Vectorization
vectorized
Vegt
vel
Velázquez
Verlag
verlet
Verlet

6020
examples/USER/misc/slater/data.after_equilibration Normal file
View File

File diff suppressed because it is too large Load Diff

6020
examples/USER/misc/slater/data.after_production_run Normal file
View File

File diff suppressed because it is too large Load Diff

60
examples/USER/misc/slater/in.slater Normal file
View File

@ -0,0 +1,60 @@
# Bulk polyelectrolyte as described in section IV of J. Chem. Phys. 125, 224107 (2006)
boundary p p p
units lj
atom_style charge
region my_sim_box block 0.0 10.0 0.0 10.0 0.0 10.0
create_box 3 my_sim_box
create_atoms 1 random 2804 100 my_sim_box
create_atoms 2 random 98 200 my_sim_box
create_atoms 3 random 98 300 my_sim_box
set type 2 charge -1.0
set type 3 charge 1.0
comm_modify mode single vel yes
mass 1 1.0
mass 2 1.0
mass 3 1.0
pair_style hybrid/overlay dpd 1.0 1.0 245455 coul/slater/long 0.929 3.0
pair_coeff * * dpd 25.0 4.5
pair_coeff * * coul/slater/long
kspace_style ewald 0.00001
dielectric 1.0
neighbor 2.0 bin
neigh_modify every 1 delay 0 check no once no
timestep 0.02
fix 2 all nve
thermo 10
thermo_style custom step spcpu temp press etotal pe ke ecoul elong evdwl
thermo_modify line one
run 100000
write_data data.after_equilibration
compute RDF_1_1 all rdf 50 1 1 cutoff 3.0
compute RDF_1_2 all rdf 50 1 2 cutoff 3.0
compute RDF_1_3 all rdf 50 1 3 cutoff 3.0
compute RDF_2_2 all rdf 50 2 2 cutoff 3.0
compute RDF_2_3 all rdf 50 2 3 cutoff 3.0
compute RDF_3_3 all rdf 50 3 3 cutoff 3.0
fix 11 all ave/time 50 1 50 c_RDF_1_1[*] file tmp_1_1.rdf mode vector
fix 12 all ave/time 50 1 50 c_RDF_1_2[*] file tmp_1_2.rdf mode vector
fix 13 all ave/time 50 1 50 c_RDF_1_3[*] file tmp_1_3.rdf mode vector
fix 14 all ave/time 50 1 50 c_RDF_2_2[*] file tmp_2_2.rdf mode vector
fix 15 all ave/time 50 1 50 c_RDF_2_3[*] file tmp_2_3.rdf mode vector
fix 16 all ave/time 50 1 50 c_RDF_3_3[*] file tmp_3_3.rdf mode vector
run 10000
write_data data.after_production_run

11225
examples/USER/misc/slater/log.lammps Normal file
View File

File diff suppressed because it is too large Load Diff

10254
examples/USER/misc/slater/tmp_1_1.rdf Normal file
View File

File diff suppressed because it is too large Load Diff

10254
examples/USER/misc/slater/tmp_1_2.rdf Normal file
View File

File diff suppressed because it is too large Load Diff

10254
examples/USER/misc/slater/tmp_1_3.rdf Normal file
View File

File diff suppressed because it is too large Load Diff

10254
examples/USER/misc/slater/tmp_2_2.rdf Normal file
View File

File diff suppressed because it is too large Load Diff

10254
examples/USER/misc/slater/tmp_2_3.rdf Normal file
View File

File diff suppressed because it is too large Load Diff

10254
examples/USER/misc/slater/tmp_3_3.rdf Normal file
View File

File diff suppressed because it is too large Load Diff

13
examples/numdiff/README.md Normal file
View File

@ -0,0 +1,13 @@
# LAMMPS FIX NUMDIFF EXAMPLE
## Numerical Difference Fix
This directory contains the ingredients to run an NVE simulation using the numerical difference fix and calculate error in forces.
Example:
```
NP=4 #number of processors
mpirun -np $NP lmp_mpi -in.numdiff
```
## Required LAMMPS packages: MOLECULE package

33
examples/numdiff/in.numdiff Normal file
View File

@ -0,0 +1,33 @@
# Numerical difference calculation of error in forces
units metal
atom_style atomic
atom_modify map yes
lattice fcc 5.358000
region box block 0 6 0 6 0 6
create_box 1 box
create_atoms 1 box
mass 1 39.903
velocity all create 10 2357 mom yes dist gaussian
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
neighbor 1 bin
timestep 0.001
fix numdiff all numdiff 200 0.0001
fix nve all nve
variable errx atom fx-f_numdiff[1]
variable erry atom fy-f_numdiff[2]
variable errz atom fz-f_numdiff[3]
write_dump all custom tmp.error f_numdiff[1] f_numdiff[2] f_numdiff[3]
dump forces all custom 200 force_error.dump v_errx v_erry v_errz
thermo 200
run 2000

16
src/GRANULAR/fix_pour.cpp
View File

@ -632,7 +632,13 @@ void FixPour::pre_exchange()
int n = atom->nlocal - 1;
atom->tag[n] = maxtag_all + m+1;
if (mode == MOLECULE) {
if (atom->molecule_flag) atom->molecule[n] = maxmol_all+1;
if (atom->molecule_flag) {
if (onemols[imol]->moleculeflag) {
atom->molecule[n] = maxmol_all + onemols[imol]->molecule[m];
} else {
atom->molecule[n] = maxmol_all+1;
}
}
if (atom->molecular == 2) {
atom->molindex[n] = 0;
atom->molatom[n] = m;
@ -666,7 +672,13 @@ void FixPour::pre_exchange()
fixshake->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
maxtag_all += natom;
if (mode == MOLECULE && atom->molecule_flag) maxmol_all++;
if (mode == MOLECULE && atom->molecule_flag) {
if (onemols[imol]->moleculeflag) {
maxmol_all += onemols[imol]->nmolecules;
} else {
maxmol_all++;
}
}
}
// warn if not successful with all insertions b/c too many attempts

16
src/MISC/fix_deposit.cpp
View File

@ -539,7 +539,13 @@ void FixDeposit::pre_exchange()
n = atom->nlocal - 1;
atom->tag[n] = maxtag_all + m+1;
if (mode == MOLECULE) {
if (atom->molecule_flag) atom->molecule[n] = maxmol_all+1;
if (atom->molecule_flag) {
if (onemols[imol]->moleculeflag) {
atom->molecule[n] = maxmol_all + onemols[imol]->molecule[m];
} else {
atom->molecule[n] = maxmol_all+1;
}
}
if (atom->molecular == 2) {
atom->molindex[n] = 0;
atom->molatom[n] = m;
@ -603,7 +609,13 @@ void FixDeposit::pre_exchange()
maxtag_all += natom;
if (maxtag_all >= MAXTAGINT)
error->all(FLERR,"New atom IDs exceed maximum allowed ID");
if (mode == MOLECULE && atom->molecule_flag) maxmol_all++;
if (mode == MOLECULE && atom->molecule_flag) {
if (onemols[imol]->moleculeflag) {
maxmol_all += onemols[imol]->nmolecules;
} else {
maxmol_all++;
}
}
if (atom->map_style) {
atom->map_init();
atom->map_set();

2
src/USER-MISC/README
View File

@ -76,6 +76,8 @@ pair_style buck/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
pair_style cosine/squared, Eugen Rozic, eugen.rozic.17 at ucl.ac.uk, 9 Aug 19
pair_style coul/diel, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11
pair_style coul/shield, Wengen Ouyang (Tel Aviv University), w.g.ouyang at gmail dot com, 30 Mar 18
pair_style coul/slater/cut, Evangelos Voyiatzis, evoyiatzis at gmail.com, 26 February 2020
pair_style coul/slater/long, Evangelos Voyiatzis, evoyiatzis at gmail.com, 26 February 2020
pair_style dipole/sf, Mario Orsi, orsimario at gmail.com, 8 Aug 11
pair_style e3b, Steven Strong (U Chicago), stevene.strong at gmail dot com, 16 Apr 19
pair_style drip, Mingjian Wen, University of Minnesota, wenxx151 at umn.edu, 17 Apr 19

187
src/USER-MISC/pair_coul_slater_cut.cpp Normal file
View File

@ -0,0 +1,187 @@
/* ----------------------------------------------------------------------
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 author: Evangelos Voyiatzis (Royal DSM)
* ------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "pair_coul_slater_cut.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairCoulSlaterCut::PairCoulSlaterCut(LAMMPS *lmp) : PairCoulCut(lmp) {}
/* ---------------------------------------------------------------------- */
void PairCoulSlaterCut::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
double rsq,r2inv,r,rinv,forcecoul,factor_coul,bracket_term;
int *ilist,*jlist,*numneigh,**firstneigh;
ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r = sqrt(rsq);
rinv = 1.0/r;
bracket_term = 1 - exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
forcecoul = qqrd2e * scale[itype][jtype] *
qtmp*q[j] * bracket_term * rinv;
fpair = factor_coul*forcecoul * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag) ecoul = factor_coul * qqrd2e *
scale[itype][jtype] * qtmp*q[j] * rinv *
(1 - (1 + r/lamda)*exp(-2*r/lamda));
if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,ecoul,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairCoulSlaterCut::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,"Illegal pair_style command");
lamda = force->numeric(FLERR,arg[0]);
cut_global = force->numeric(FLERR,arg[1]);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCoulSlaterCut::write_restart_settings(FILE *fp)
{
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&lamda,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCoulSlaterCut::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_global,sizeof(double),1,fp);
fread(&lamda,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&lamda,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
double PairCoulSlaterCut::single(int i, int j, int /*itype*/, int /*jtype*/,
double rsq, double factor_coul, double /*factor_lj*/,
double &fforce)
{
double r2inv,r,rinv,forcecoul,phicoul,bracket_term;
r2inv = 1.0/rsq;
r = sqrt(rsq);
rinv = 1.0/r;
bracket_term = 1 - exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
forcecoul = force->qqrd2e * atom->q[i]*atom->q[j] *
bracket_term * rinv;
fforce = factor_coul*forcecoul * r2inv;
phicoul = force->qqrd2e * atom->q[i]*atom->q[j] * rinv * (1 - (1 + r/lamda)*exp(-2*r/lamda));
return factor_coul*phicoul;
}

53
src/USER-MISC/pair_coul_slater_cut.h Normal file
View File

@ -0,0 +1,53 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(coul/slater/cut,PairCoulSlaterCut)
#else
#ifndef LMP_PAIR_COUL_SLATER_CUT_H
#define LMP_PAIR_COUL_SLATER_CUT_H
#include "pair_coul_cut.h"
namespace LAMMPS_NS {
class PairCoulSlaterCut : public PairCoulCut {
public:
PairCoulSlaterCut(class LAMMPS *);
virtual void compute(int, int);
void settings(int, char **);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
protected:
double lamda;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/

422
src/USER-MISC/pair_coul_slater_long.cpp Normal file
View File

@ -0,0 +1,422 @@
/* ----------------------------------------------------------------------
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 author: Evangelos Voyiatzis (Royal DSM)
* ------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "pair_coul_slater_long.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "kspace.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "update.h"
#include "integrate.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
/* ---------------------------------------------------------------------- */
PairCoulSlaterLong::PairCoulSlaterLong(LAMMPS *lmp) : Pair(lmp)
{
ewaldflag = pppmflag = 1;
//ftable = NULL;
qdist = 0.0;
}
/* ---------------------------------------------------------------------- */
PairCoulSlaterLong::~PairCoulSlaterLong()
{
if (!copymode) {
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(scale);
}
//if (ftable) free_tables();
}
}
/* ---------------------------------------------------------------------- */
void PairCoulSlaterLong::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
// double fraction,table;
double r,r2inv,forcecoul,factor_coul;
double grij,expm2,prefactor,t,erfc;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq;
double slater_term;
// int itable;
ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cut_coulsq) {
r2inv = 1.0/rsq;
// if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
slater_term = exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - slater_term);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
/*
} else {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
if (factor_coul < 1.0) {
table = ctable[itable] + fraction*dctable[itable];
prefactor = scale[itype][jtype] * qtmp*q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
}
}
*/
fpair = forcecoul * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag) {
// if (!ncoultablebits || rsq <= tabinnersq)
ecoul = prefactor*(erfc - (1 + r/lamda)*exp(-2*r/lamda));
/*
else {
table = etable[itable] + fraction*detable[itable];
ecoul = scale[itype][jtype] * qtmp*q[j] * table;
}
*/
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,ecoul,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairCoulSlaterLong::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(scale,n+1,n+1,"pair:scale");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairCoulSlaterLong::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,"Illegal pair_style command");
lamda = force->numeric(FLERR,arg[0]);
cut_coul = force->numeric(FLERR,arg[1]);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairCoulSlaterLong::coeff(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
scale[i][j] = 1.0;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairCoulSlaterLong::init_style()
{
if (!atom->q_flag)
error->all(FLERR,"Pair style coul/slater/long requires atom attribute q");
neighbor->request(this,instance_me);
cut_coulsq = cut_coul * cut_coul;
// insure use of KSpace long-range solver, set g_ewald
if (force->kspace == NULL)
error->all(FLERR,"Pair style requires a KSpace style");
g_ewald = force->kspace->g_ewald;
// setup force tables
// if (ncoultablebits) init_tables(cut_coul,NULL);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairCoulSlaterLong::init_one(int i, int j)
{
scale[j][i] = scale[i][j];
return cut_coul+2.0*qdist;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCoulSlaterLong::write_restart(FILE *fp)
{
write_restart_settings(fp);
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j])
fwrite(&scale[i][j],sizeof(double),1,fp);
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCoulSlaterLong::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) fread(&scale[i][j],sizeof(double),1,fp);
MPI_Bcast(&scale[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCoulSlaterLong::write_restart_settings(FILE *fp)
{
fwrite(&cut_coul,sizeof(double),1,fp);
fwrite(&lamda,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
//fwrite(&ncoultablebits,sizeof(int),1,fp);
//fwrite(&tabinner,sizeof(double),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCoulSlaterLong::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_coul,sizeof(double),1,fp);
fread(&lamda,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
//fread(&ncoultablebits,sizeof(int),1,fp);
//fread(&tabinner,sizeof(double),1,fp);
}
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
MPI_Bcast(&lamda,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
//MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
//MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
}
/* ---------------------------------------------------------------------- */
double PairCoulSlaterLong::single(int i, int j, int /*itype*/, int /*jtype*/,
double rsq,
double factor_coul, double /*factor_lj*/,
double &fforce)
{
double r2inv,r,grij,expm2,t,erfc,prefactor;
double slater_term;
// double fraction,table;
double forcecoul,phicoul;
// int itable;
r2inv = 1.0/rsq;
// if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
slater_term = exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - slater_term);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
/*
} else {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = atom->q[i]*atom->q[j] * table;
if (factor_coul < 1.0) {
table = ctable[itable] + fraction*dctable[itable];
prefactor = atom->q[i]*atom->q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
}
}
*/
fforce = forcecoul * r2inv;
// if (!ncoultablebits || rsq <= tabinnersq)
phicoul = prefactor*(erfc - (1 + r/lamda)*exp(-2*r/lamda));
/*
else {
table = etable[itable] + fraction*detable[itable];
phicoul = atom->q[i]*atom->q[j] * table;
}
*/
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
return phicoul;
}
/* ---------------------------------------------------------------------- */
void *PairCoulSlaterLong::extract(const char *str, int &dim)
{
if (strcmp(str,"cut_coul") == 0) {
dim = 0;
return (void *) &cut_coul;
}
if (strcmp(str,"lamda") == 0) {
dim = 0;
return (void *) &lamda;
}
if (strcmp(str,"scale") == 0) {
dim = 2;
return (void *) scale;
}
return NULL;
}

78
src/USER-MISC/pair_coul_slater_long.h Normal file
View File

@ -0,0 +1,78 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(coul/slater/long,PairCoulSlaterLong)
#else
#ifndef LMP_PAIR_COUL_SLATER_LONG_H
#define LMP_PAIR_COUL_SLATER_LONG_H
#include "pair.h"
namespace LAMMPS_NS {
class PairCoulSlaterLong : public Pair {
public:
PairCoulSlaterLong(class LAMMPS *);
~PairCoulSlaterLong();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
virtual void init_style();
virtual double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual double single(int, int, int, int, double, double, double, double &);
virtual void *extract(const char *, int &);
protected:
double cut_coul,cut_coulsq,qdist;
double lamda;
//double *cut_respa;
double g_ewald;
double **scale;
virtual void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair style coul/slater/long requires atom attribute q
The atom style defined does not have this attribute.
E: Pair style requires a KSpace style
No kspace style is defined.
*/

15
src/create_atoms.cpp
View File

@ -490,7 +490,13 @@ void CreateAtoms::command(int narg, char **arg)
for (int i = 0; i < molcreate; i++) {
if (tag) offset = tag[ilocal]-1;
for (int m = 0; m < natoms; m++) {
if (molecule_flag) molecule[ilocal] = moloffset + i+1;
if (molecule_flag) {
if (onemol->moleculeflag) {
molecule[ilocal] = moloffset + onemol->molecule[m];
} else {
molecule[ilocal] = moloffset + 1;
}
}
if (molecular == 2) {
atom->molindex[ilocal] = 0;
atom->molatom[ilocal] = m;
@ -524,6 +530,13 @@ void CreateAtoms::command(int narg, char **arg)
}
ilocal++;
}
if (molecule_flag) {
if (onemol->moleculeflag) {
moloffset += onemol->nmolecules;
} else {
moloffset++;
}
}
}
// perform irregular comm to migrate atoms to new owning procs

352
src/fix_numdiff.cpp Normal file
View File

@ -0,0 +1,352 @@
/* ----------------------------------------------------------------------
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 author: Charles Sievers (UC Davis)
------------------------------------------------------------------------- */
#include "fix_numdiff.h"
#include <mpi.h>
#include <cstring>
#include "atom.h"
#include "domain.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "respa.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixNumDiff::FixNumDiff(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), numdiff_forces(NULL), temp_f(NULL),
id_pe(NULL)
{
if (narg < 5) error->all(FLERR,"Illegal fix numdiff command");
peratom_flag = 1;
peratom_freq = nevery;
size_peratom_cols = 3;
respa_level_support = 1;
nevery = force->inumeric(FLERR,arg[3]);
delta = force->numeric(FLERR,arg[4]);
if (nevery <= 0 || delta <= 0.0)
error->all(FLERR,"Illegal fix numdiff command");
int n = strlen(id) + 6;
id_pe = new char[n];
strcpy(id_pe,id);
strcat(id_pe,"_pe");
char **newarg = new char*[3];
newarg[0] = id_pe;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pe";
modify->add_compute(3,newarg);
delete [] newarg;
maxatom = 0;
numdiff_forces = NULL;
temp_x = NULL;
temp_f = NULL;
if (atom->map_style == 0)
error->all(FLERR,"Fix numdiff requires an atom map, see atom_modify");
// perform initial allocation of atom-based arrays
// zero numdiff_forces since dump may access it on timestep 0
// zero numdiff_forces since a variable may access it before first run
reallocate();
force_clear(numdiff_forces);
}
/* ---------------------------------------------------------------------- */
FixNumDiff::~FixNumDiff()
{
memory->destroy(numdiff_forces);
memory->destroy(temp_x);
memory->destroy(temp_f);
modify->delete_compute(id_pe);
delete [] id_pe;
}
/* ---------------------------------------------------------------------- */
int FixNumDiff::setmask()
{
datamask_read = datamask_modify = 0;
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixNumDiff::init()
{
// require consecutive atom IDs
if (!atom->tag_enable || !atom->tag_consecutive())
error->all(FLERR,"Fix numdiff requires consecutive atom IDs");
// check for PE compute
int icompute = modify->find_compute(id_pe);
if (icompute < 0) error->all(FLERR,"Compute ID for fix numdiff does not exist");
pe = modify->compute[icompute];
if (force->pair && force->pair->compute_flag) pair_compute_flag = 1;
else pair_compute_flag = 0;
if (force->kspace && force->kspace->compute_flag) kspace_compute_flag = 1;
else kspace_compute_flag = 0;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
void FixNumDiff::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
void FixNumDiff::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixNumDiff::post_force(int vflag)
{
if (update->ntimestep % nevery) return;
calculate_forces();
}
/* ---------------------------------------------------------------------- */
void FixNumDiff::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixNumDiff::min_post_force(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
compute finite difference forces
------------------------------------------------------------------------- */
void FixNumDiff::calculate_forces()
{
int i,j,ilocal;
double energy;
// grow arrays if necessary
if (atom->nlocal + atom->nghost > maxatom) reallocate();
// store copy of current forces for owned and ghost atoms
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
for (i = 0; i < nall; i++)
for (j = 0; j < 3; j++) {
temp_x[i][j] = x[i][j];
temp_f[i][j] = f[i][j];
}
// initialize numerical forces to zero
force_clear(numdiff_forces);
// loop over all atoms in system
// compute a finite difference force in each dimension
int flag,allflag;
double denominator = 0.5 / delta;
int *mask = atom->mask;
int ntotal = static_cast<tagint> (atom->natoms);
int dimension = domain->dimension;
for (tagint m = 1; m <= ntotal; m++) {
ilocal = atom->map(m);
flag = 0;
if ((ilocal >= 0 && ilocal < nlocal) && (mask[ilocal] & groupbit)) flag = 1;
MPI_Allreduce(&flag,&allflag,1,MPI_INT,MPI_SUM,world);
if (!allflag) continue;
for (int idim = 0; idim < dimension; idim++) {
displace_atoms(ilocal,idim,1);
energy = update_energy();
if (ilocal >= 0 && ilocal < nlocal)
numdiff_forces[ilocal][idim] -= energy;
displace_atoms(ilocal,idim,-2);
energy = update_energy();
if (ilocal >= 0 && ilocal < nlocal) {
numdiff_forces[ilocal][idim] += energy;
numdiff_forces[ilocal][idim] *= denominator;
}
restore_atoms(ilocal,idim);
}
}
// restore original forces for owned and ghost atoms
for (i = 0; i < nall; i++)
for (j = 0; j < 3; j++) {
f[i][j] = temp_f[i][j];
}
}
/* ----------------------------------------------------------------------
displace position of all owned and ghost copies of ilocal
---------------------------------------------------------------------- */
void FixNumDiff::displace_atoms(int ilocal, int idim, int magnitude)
{
if (ilocal < 0) return;
double **x = atom->x;
int *sametag = atom->sametag;
int j = ilocal;
x[ilocal][idim] += delta*magnitude;
while (sametag[j] >= 0) {
j = sametag[j];
x[j][idim] += delta*magnitude;
}
}
/* ----------------------------------------------------------------------
restore position of all owned and ghost copies of ilocal
---------------------------------------------------------------------- */
void FixNumDiff::restore_atoms(int ilocal, int idim)
{
if (ilocal < 0) return;
double **x = atom->x;
int *sametag = atom->sametag;
int j = ilocal;
x[ilocal][idim] = temp_x[ilocal][idim];
while (sametag[j] >= 0) {
j = sametag[j];
x[j][idim] = temp_x[j][idim];
}
}
/* ----------------------------------------------------------------------
evaluate potential energy and forces
same logic as in Verlet
------------------------------------------------------------------------- */
double FixNumDiff::update_energy()
{
force_clear(atom->f);
int eflag = 1;
if (pair_compute_flag) force->pair->compute(eflag,0);
if (atom->molecular) {
if (force->bond) force->bond->compute(eflag,0);
if (force->angle) force->angle->compute(eflag,0);
if (force->dihedral) force->dihedral->compute(eflag,0);
if (force->improper) force->improper->compute(eflag,0);
}
if (kspace_compute_flag) force->kspace->compute(eflag,0);
double energy = pe->compute_scalar();
return energy;
}
/* ----------------------------------------------------------------------
clear forces needed
------------------------------------------------------------------------- */
void FixNumDiff::force_clear(double **forces)
{
size_t nbytes = sizeof(double) * atom->nlocal;
if (force->newton) nbytes += sizeof(double) * atom->nghost;
if (nbytes) memset(&forces[0][0],0,3*nbytes);
}
/* ----------------------------------------------------------------------
reallocated local per-atoms arrays
------------------------------------------------------------------------- */
void FixNumDiff::reallocate()
{
memory->destroy(numdiff_forces);
memory->destroy(temp_x);
memory->destroy(temp_f);
maxatom = atom->nmax;
memory->create(numdiff_forces,maxatom,3,"numdiff:numdiff_force");
memory->create(temp_x,maxatom,3,"numdiff:temp_x");
memory->create(temp_f,maxatom,3,"numdiff:temp_f");
array_atom = numdiff_forces;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixNumDiff::memory_usage()
{
bigint bytes = 0.0;
bytes += 3 * maxatom*3 * sizeof(double);
return bytes;
}

90
src/fix_numdiff.h Normal file
View File

@ -0,0 +1,90 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(numdiff,FixNumDiff)
#else
#ifndef LMP_FIX_NUMDIFF_H
#define LMP_FIX_NUMDIFF_H
#include "fix.h"
namespace LAMMPS_NS {
class FixNumDiff : public Fix {
public:
FixNumDiff(class LAMMPS *, int, char **);
~FixNumDiff();
int setmask();
void init();
void setup(int);
void min_setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
double memory_usage();
private:
double delta;
int maxatom;
int ilevel_respa;
int pair_compute_flag; // 0 if pair->compute is skipped
int kspace_compute_flag; // 0 if kspace->compute is skipped
char *id_pe;
class Compute *pe;
double **numdiff_forces; // finite diff forces
double **temp_x; // original coords
double **temp_f; // original forces
void calculate_forces();
void displace_atoms(int, int, int);
void restore_atoms(int, int);
double update_energy();
void force_clear(double **);
void reallocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix numdiff requires an atom map, see atom_modify
Self-explanatory. Efficient loop over all atoms for numerical difference
requires an atom map.
E: Fix numdiff requires consecutive atom IDs
Self-explanatory. Efficient loop over all atoms for numerical difference
requires consecutive atom IDs.
E: Compute ID for fix numdiff does not exist
Self-explanatory.
*/

101
src/molecule.cpp
View File

@ -37,7 +37,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
Pointers(lmp), id(NULL), x(NULL), type(NULL), q(NULL), radius(NULL),
Pointers(lmp), id(NULL), x(NULL), type(NULL), molecule(NULL), q(NULL), radius(NULL),
rmass(NULL), num_bond(NULL), bond_type(NULL), bond_atom(NULL),
num_angle(NULL), angle_type(NULL), angle_atom1(NULL), angle_atom2(NULL),
angle_atom3(NULL), num_dihedral(NULL), dihedral_type(NULL), dihedral_atom1(NULL),
@ -46,7 +46,7 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
improper_atom3(NULL), improper_atom4(NULL), nspecial(NULL), special(NULL),
shake_flag(NULL), shake_atom(NULL), shake_type(NULL), avec_body(NULL), ibodyparams(NULL),
dbodyparams(NULL), dx(NULL), dxcom(NULL), dxbody(NULL), quat_external(NULL),
fp(NULL), count(NULL)
fp(NULL), count(NULL), fragmentmask(NULL)
{
me = comm->me;
@ -146,19 +146,19 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
if (me == 0) {
if (screen)
fprintf(screen,"Read molecule %s:\n"
fprintf(screen,"Read molecule template %s:\n %d molecules\n"
" %d atoms with max type %d\n %d bonds with max type %d\n"
" %d angles with max type %d\n %d dihedrals with max type %d\n"
" %d impropers with max type %d\n",
id,natoms,ntypes,
id,nmolecules,natoms,ntypes,
nbonds,nbondtypes,nangles,nangletypes,
ndihedrals,ndihedraltypes,nimpropers,nimpropertypes);
if (logfile)
fprintf(logfile,"Read molecule %s:\n"
fprintf(logfile,"Read molecule template %s:\n %d molecules\n"
" %d atoms with max type %d\n %d bonds with max type %d\n"
" %d angles with max type %d\n %d dihedrals with max type %d\n"
" %d impropers with max type %d\n",
id,natoms,ntypes,
id,nmolecules,natoms,ntypes,
nbonds,nbondtypes,nangles,nangletypes,
ndihedrals,ndihedraltypes,nimpropers,nimpropertypes);
}
@ -447,6 +447,9 @@ void Molecule::read(int flag)
} else if (strstr(line,"impropers")) {
nmatch = sscanf(line,"%d",&nimpropers);
nwant = 1;
} else if (strstr(line,"fragments")) {
nmatch = sscanf(line,"%d",&nfragments);
nwant = 1;
} else if (strstr(line,"mass")) {
massflag = 1;
nmatch = sscanf(line,"%lg",&masstotal);
@ -519,6 +522,14 @@ void Molecule::read(int flag)
typeflag = 1;
if (flag) types(line);
else skip_lines(natoms,line);
} else if (strcmp(keyword,"Molecules") == 0) {
moleculeflag = 1;
if (flag) molecules(line);
else skip_lines(natoms,line);
} else if (strcmp(keyword,"Fragments") == 0) {
fragmentflag = 1;
if (flag) fragments(line);
else skip_lines(nfragments,line);
} else if (strcmp(keyword,"Charges") == 0) {
qflag = 1;
if (flag) charges(line);
@ -694,6 +705,58 @@ void Molecule::types(char *line)
ntypes = MAX(ntypes,type[i]);
}
/* ----------------------------------------------------------------------
read molecules from file
set nmolecules = max of any molecule type
------------------------------------------------------------------------- */
void Molecule::molecules(char *line)
{
int tmp;
for (int i = 0; i < natoms; i++) {
readline(line);
if (2 != sscanf(line,"%d %d",&tmp,&molecule[i]))
error->all(FLERR,"Invalid Molecules section in molecule file");
// molecule[i] += moffset; // placeholder for possible molecule offset
}
for (int i = 0; i < natoms; i++)
if (molecule[i] <= 0)
error->all(FLERR,"Invalid molecule ID in molecule file");
for (int i = 0; i < natoms; i++)
nmolecules = MAX(nmolecules,molecule[i]);
}
/* ----------------------------------------------------------------------
read fragments from file
------------------------------------------------------------------------- */
void Molecule::fragments(char *line)
{
int n,m,atomID,nwords;
char **words = new char*[natoms+1];
for (int i = 0; i < nfragments; i++) {
readline(line);
nwords = parse(line,words,natoms+1);
if (nwords > natoms+1)
error->all(FLERR,"Invalid atom ID in Fragments section of molecule file");
n = strlen(words[0]) + 1;
fragmentnames[i] = new char[n];
strcpy(fragmentnames[i],words[0]);
for (m = 1; m < nwords; m++) {
atomID = atoi(words[m]);
if (atomID <= 0 || atomID > natoms)
error->all(FLERR,"Invalid atom ID in Fragments section of molecule file");
fragmentmask[i][atomID-1] = 1;
}
}
delete [] words;
}
/* ----------------------------------------------------------------------
read charges from file
------------------------------------------------------------------------- */
@ -1359,6 +1422,17 @@ void Molecule::body(int flag, int pflag, char *line)
}
}
/* ----------------------------------------------------------------------
return fragment index if name matches existing fragment, -1 if no such fragment
------------------------------------------------------------------------- */
int Molecule::findfragment(const char *name)
{
for (int i = 0; i < nfragments; i++)
if (fragmentnames[i] && strcmp(name,fragmentnames[i]) == 0) return i;
return -1;
}
/* ----------------------------------------------------------------------
error check molecule attributes and topology against system settings
flag = 0, just check this molecule
@ -1432,13 +1506,14 @@ void Molecule::initialize()
natoms = 0;
nbonds = nangles = ndihedrals = nimpropers = 0;
ntypes = 0;
nmolecules = 1;
nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
nibody = ndbody = 0;
bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
maxspecial = 0;
xflag = typeflag = qflag = radiusflag = rmassflag = 0;
xflag = typeflag = moleculeflag = fragmentflag = qflag = radiusflag = rmassflag = 0;
bondflag = angleflag = dihedralflag = improperflag = 0;
nspecialflag = specialflag = 0;
shakeflag = shakeflagflag = shakeatomflag = shaketypeflag = 0;
@ -1493,6 +1568,11 @@ void Molecule::allocate()
{
if (xflag) memory->create(x,natoms,3,"molecule:x");
if (typeflag) memory->create(type,natoms,"molecule:type");
if (moleculeflag) memory->create(molecule,natoms,"molecule:molecule");
if (fragmentflag) fragmentnames = new char*[nfragments];
if (fragmentflag) memory->create(fragmentmask,nfragments,natoms,"molecule:fragmentmask");
for (int i = 0; i < nfragments; i++)
for (int j = 0; j < natoms; j++) fragmentmask[i][j] = 0;
if (qflag) memory->create(q,natoms,"molecule:q");
if (radiusflag) memory->create(radius,natoms,"molecule:radius");
if (rmassflag) memory->create(rmass,natoms,"molecule:rmass");
@ -1580,10 +1660,17 @@ void Molecule::deallocate()
{
memory->destroy(x);
memory->destroy(type);
memory->destroy(molecule);
memory->destroy(q);
memory->destroy(radius);
memory->destroy(rmass);
memory->destroy(fragmentmask);
if (fragmentflag) {
for (int i = 0; i < nfragments; i++) delete [] fragmentnames[i];
delete [] fragmentnames;
}
memory->destroy(num_bond);
memory->destroy(bond_type);
memory->destroy(bond_atom);

24
src/molecule.h
View File

@ -30,10 +30,14 @@ class Molecule : protected Pointers {
int natoms;
int nbonds,nangles,ndihedrals,nimpropers;
int ntypes;
int ntypes,nmolecules,nfragments;
int nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;
int nibody,ndbody;
// fragment info
char **fragmentnames;
int **fragmentmask; // nfragments by natoms
// max bond,angle,etc per atom
int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
@ -41,7 +45,7 @@ class Molecule : protected Pointers {
// 1 if attribute defined in file, 0 if not
int xflag,typeflag,qflag,radiusflag,rmassflag;
int xflag,typeflag,moleculeflag,fragmentflag,qflag,radiusflag,rmassflag;
int bondflag,angleflag,dihedralflag,improperflag;
int nspecialflag,specialflag;
int shakeflag,shakeflagflag,shakeatomflag,shaketypeflag;
@ -59,6 +63,7 @@ class Molecule : protected Pointers {
double **x; // displacement of each atom from origin
int *type; // type of each atom
tagint *molecule; // molecule of each atom
double *q; // charge on each atom
double *radius; // radius of each atom
double *rmass; // mass of each atom
@ -118,6 +123,7 @@ class Molecule : protected Pointers {
void compute_mass();
void compute_com();
void compute_inertia();
int findfragment(const char *);
void check_attributes(int);
private:
@ -131,6 +137,8 @@ class Molecule : protected Pointers {
void read(int);
void coords(char *);
void types(char *);
void molecules(char *);
void fragments(char *);
void charges(char *);
void diameters(char *);
void masses(char *);
@ -363,6 +371,18 @@ E: Invalid improper type in impropers section of molecule file
Self-explanatory.
E: Invalid molecule ID in molecule file
Molecule ID must be a non-zero positive integer.
E: Invalid Molecules section in molecule file
Self-explanatory.
E: Invalid atom ID in Fragments section of molecule file
Self-explanatory.
E: Invalid Special Bond Counts section in molecule file
Self-explanatory.