ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #1503 from etomica/master

Browse Source 
Implement HMA compute in LAMMPS
This commit is contained in:
Axel Kohlmeyer
2019-08-21 00:12:16 -04:00
committed by GitHub
parent 42b282fe36 91a80c3c21
commit cfa9179b28
16 changed files with 21209 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/compute.txt
View File

@ -217,6 +217,7 @@ compute"_Commands_compute.html doc page are followed by one or more of
"heat/flux"_compute_heat_flux.html - heat flux through a group of atoms
"heat/flux/tally"_compute_tally.html -
"hexorder/atom"_compute_hexorder_atom.html - bond orientational order parameter q6
"hma"_compute_hma.html - harmonically mapped averaging for atomic crystals
"improper"_compute_improper.html - energy of each improper sub-style
"improper/local"_compute_improper_local.html - angle of each improper
"inertia/chunk"_compute_inertia_chunk.html - inertia tensor for each chunk

184
doc/src/compute_hma.txt Normal file
View File

@ -0,0 +1,184 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
compute hma command :h3
[Syntax:]
compute ID group-ID hma temp-ID keyword ... :pre
ID, group-ID are documented in "compute"_compute.html command :l
hma = style name of this compute command :l
temp-ID = ID of fix that specifies the set temperature during canonical simulation :l
keyword = {anharmonic} {u} {p Pharm} {cv} :l
{anharmonic} = compute will return anharmonic property values
{u} = compute will return potential energy
{p} = compute will return pressure. the following keyword must be the difference between the harmonic pressure and lattice pressure as described below
{cv} = compute will return the heat capacity :pre
:ule
[Examples:]
compute 2 all hma 1 u
compute 2 all hma 1 anharmonic u p 0.9
compute 2 all hma 1 u cv :pre
[Description:]
Define a computation that calculates the properties of a solid (potential
energy, pressure or heat capacity), using the harmonically-mapped averaging
(HMA) method.
This command yields much higher precision than the equivalent compute commands
("compute pe"_compute_pe.html, "compute pressure"_compute_pressure.html, etc.)
commands during a canonical simulation of an atomic crystal. Specifically,
near melting HMA can yield averages of a given precision an order of magnitude
faster than conventional methods, and this only improves as the temperatures is
lowered. This is particularly important for evaluating the free energy by
thermodynamic integration, where the low-temperature contributions are the
greatest source of statistical uncertainty. Moreover, HMA has other
advantages, including smaller potential-truncation effects, finite-size
effects, smaller timestep inaccuracy, faster equilibration and shorter
decorrelation time.
HMA should not be used if atoms are expected to diffuse. It is also
restricted to simulations in the NVT ensemble. While this compute may be
used with any potential in LAMMPS, it will provide inaccurate results
for potentials that do not go to 0 at the truncation distance;
"pair_lj_smooth_linear"_pair_lj_smooth_linear.html and Ewald summation should
work fine, while "pair_lj"_pair_lj.html will perform poorly unless
the potential is shifted (via "pair_modify"_pair_modify.html shift) or the cutoff is large. Furthermore, computation of the heat capacity with
this compute is restricted to those that implement the single_hessian method
in Pair. Implementing single_hessian in additional pair styles is simple.
Please contact Andrew Schultz (ajs42 at buffalo.edu) and David Kofke (kofke at
buffalo.edu) if your desired pair style does not have this method. This is
the list of pair styles that currently implement pair_hessian:
"lj_smooth_linear"_pair_lj_smooth_linear.html :l
:ule
In this method, the analytically known harmonic behavior of a crystal is removed from the traditional ensemble
averages, which leads to an accurate and precise measurement of the anharmonic contributions without contamination
by noise produced by the already-known harmonic behavior.
A detailed description of this method can be found in ("Moustafa"_#hma-Moustafa). The potential energy is computed by the formula:
\begin\{equation\}
\left< U\right>_\{HMA\} = \frac\{d\}\{2\} (N-1) k_B T + \left< U + \frac\{1\}\{2\} F\bullet\Delta r \right>
\end\{equation\}
where \(N\) is the number of atoms in the system, \(k_B\) is Boltzmann's
constant, \(T\) is the temperature, \(d\) is the
dimensionality of the system (2 or 3 for 2d/3d), \(F\bullet\Delta r\) is the sum of dot products of the
atomic force vectors and displacement (from lattice sites) vectors, and \(U\) is the sum of
pair, bond, angle, dihedral, improper, kspace (long-range), and fix energies.
The pressure is computed by the formula:
\begin\{equation\}
\left< P\right>_\{HMA\} = \Delta \hat P + \left< P_\{vir\} + \frac\{\beta \Delta \hat P - \rho\}\{d(N-1)\} F\bullet\Delta r \right>
\end\{equation\}
where \(\rho\) is the number density of the system, \(\Delta \hat P\) is the
difference between the harmonic and lattice pressure, \(P_\{vir\}\) is
the virial pressure computed as the sum of pair, bond, angle, dihedral,
improper, kspace (long-range), and fix contributions to the force on each
atom, and \(k_B=1/k_B T\). Although the method will work for any value of \(\Delta \hat P\)
specified (use pressure "units"_units.html), the precision of the resultant
pressure is sensitive to \(\Delta \hat P\); the precision tends to be
best when \(\Delta \hat P\) is the actual the difference between the lattice
pressure and harmonic pressure.
\begin\{equation\}
\left<C_V \right>_\{HMA\} = \frac\{d\}\{2\} (N-1) k_B + \frac\{1\}\{k_B T^2\} \left( \left<
U_\{HMA\}^2 \right> - \left<U_\{HMA\}\right>^2 \right) + \frac\{1\}\{4 T\}
\left< F\bullet\Delta r + \Delta r \bullet \Phi \bullet \Delta r \right>
\end\{equation\}
where \(\Phi\) is the Hessian matrix. The compute hma command
computes the full expression for \(C_V\) except for the
\(\left<U_\{HMA\}^2\right>^2\) in the variance term, which can be obtained by
passing the {u} keyword; you must add this extra contribution to the \(C_V\)
value reported by this compute. The variance term can cause significant
round-off error when computing \(C_V\). To address this, the {anharmonic}
keyword can be passed and/or the output format can be specified with more
digits.
thermo_modify format float '%22.15e' :pre
The {anharmonic} keyword will instruct the compute to return anharmonic
properties rather than the full properties, which include lattice, harmonic
and anharmonic contributions.
When using this keyword, the compute must be first active (it must be included
via a "thermo_style custom"_thermo_style.html command) while the atoms are
still at their lattice sites (before equilibration).
The temp-ID specified with compute hma command should be same as the fix-ID of Nose-Hoover ("fix nvt"_fix_nh.html) or
Berendsen ("fix temp/berendsen"_fix_temp_berendsen.html) thermostat used for the simulation. While using this command, Langevin thermostat
("fix langevin"_fix_langevin.html)
should be avoided as its extra forces interfere with the HMA implementation.
NOTE: Compute hma command should be used right after the energy minimization, when the atoms are at their lattice sites.
The simulation should not be started before this command has been used in the input script.
The following example illustrates the placement of this command in the input script:
min_style cg
minimize 1e-35 1e-15 50000 500000
compute 1 all hma thermostatid u
fix thermostatid all nvt temp 600.0 600.0 100.0 :pre
NOTE: Compute hma should be used when the atoms of the solid do not diffuse. Diffusion will reduce the precision in the potential energy computation.
NOTE: The "fix_modify energy yes"_fix_modify.html command must also be specified if a fix is to contribute potential energy to this command.
An example input script that uses this compute is included in
examples/USER/hma/ along with corresponding LAMMPS output showing that the HMA
properties fluctuate less than the corresponding conventional properties.
[Output info:]
This compute calculates a global vector that includes the n properties
requested as arguments to the command (the potential energy, pressure and/or heat
capacity). The elements of the vector can be accessed by indices 1-n by any
command that uses global vector values as input. See the "Howto
output"_Howto_output.html doc page for an overview of LAMMPS output options.
The vector values calculated by this compute are "extensive". The
scalar value will be in energy "units"_units.html.
[Restrictions:]
This compute is part of the USER-MISC package. It is enabled only
if LAMMPS was built with that package. See the "Build
package"_Build_package.html doc page for more info.
Usage restricted to canonical (NVT) ensemble simulation only.
[Related commands:]
"compute pe"_compute_pe.html, "compute pressure"_compute_pressure.html
"dynamical matrix"_dynamical_matrix.html provides a finite difference
formulation of the hessian provided by Pair's single_hessian, which is used by
this compute.
[Default:] none
:line
:link(hma-Moustafa)
[(Moustafa)] Sabry G. Moustafa, Andrew J. Schultz, and David A. Kofke, {Very fast averaging of thermal properties of crystals by molecular simulation},
"Phys. Rev. E \[92\], 043303 (2015)"_https://link.aps.org/doi/10.1103/PhysRevE.92.043303

1
doc/src/computes.txt
View File

@ -47,6 +47,7 @@ Computes :h1
compute_gyration_shape
compute_heat_flux
compute_hexorder_atom
compute_hma
compute_improper
compute_improper_local
compute_inertia_chunk

3
doc/src/dynamical_matrix.txt
View File

@ -47,6 +47,9 @@ package"_Build_package.html doc page for more info.
"fix phonon"_fix_phonon.html
"compute hma"_compute_hma.html uses an analytic formulation of the hessian
provided by Pair's single_hessian.
[Default:]
The default settings are file = "dynmat.dyn", binary = no

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

@ -45,6 +45,7 @@ Aidan
aij
airebo
Aj
ajs
ajaramil
akohlmey
Aktulga
@ -99,6 +100,7 @@ antisymmetry
anton
Antonelli
api
Apoorva
Appl
Apu
arccos
@ -523,6 +525,7 @@ Dcut
de
dE
De
decorrelation
debye
Debye
Decius
@ -1075,6 +1078,7 @@ Hilger
histo
histogrammed
histogramming
hma
hmaktulga
hoc
Hochbruck
@ -1331,6 +1335,8 @@ kmax
Kmax
Knizhnik
knl
Kofke
kofke
Kohlmeyer
Kohn
kokkos
@ -1718,6 +1724,7 @@ Morteza
Mosayebi
Moseler
Moskalev
Moustafa
mov
mpi
MPI
@ -2217,6 +2224,7 @@ PTM
ptr
pu
purdue
Purohit
pushstore
pvar
pw
@ -2434,6 +2442,7 @@ Ryckaert
Rycroft
Rydbergs
Rz
Sabry
saddlebrown
Sadigh
saed

22
examples/USER/misc/hma/README Normal file
View File

@ -0,0 +1,22 @@
The example input script sets up a simple FCC crystal using the Lennard-Jones
potential. The script sets up the HMA compute to calculate the energy, pressure
and heat capacity. The output columns are:
1: timestep
2: measured temperature
3: potential energy
4: HMA potential energy
5: pressure
6: HMA pressure
7: HMA heat capacity contribution
Averages of the potential energy (#3 and #4) agree although #4 (HMA) is more precise.
Averages of the pressure (#5 and #6) agree once the ideal gas
contribution is included; #6 (HMA) is more precise.
The heat capacity can be computed from colume #3 (convential) as
Cv = Var(#3)/(k T^2)
With HMA, the heat capacity can be computed from column #4 and #7 as
Cv = #7 + Var(#4)/(k T^2)

36
examples/USER/misc/hma/hma.in Normal file
View File

@ -0,0 +1,36 @@
# Harmonically mapped average example
units lj
dimension 3
boundary p p p
atom_style atomic
atom_modify map array
# ---------- Create Atoms ----------------------------
lattice fcc 1.0
region box block 0 4 0 4 0 4 units lattice
create_box 1 box
lattice fcc 1.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 1 region box
# ---------- Define Interatomic Potential ---------------------
pair_style lj/smooth/linear 3
pair_coeff * * 1.0 1.0
mass 1 1.0
atom_modify sort 0 1
velocity all create 0.1 45678 dist gaussian
compute u all pe
compute p all pressure NULL pair
compute hma all hma settemp u p 9.579586686264458 cv
timestep 0.005
fix settemp all nvt temp 1.0 1.0 0.5
thermo_style custom elapsed temp c_u c_hma[1] c_p c_hma[2] c_hma[3]
thermo_modify format float '%22.15e'
thermo 500
run 20000
thermo 20
run 200000

10159
examples/USER/misc/hma/log.6Nov18.hma.g++.1 Normal file
View File

File diff suppressed because it is too large Load Diff

10159
examples/USER/misc/hma/log.6Nov18.hma.g++.4 Normal file
View File

File diff suppressed because it is too large Load Diff

1
src/USER-MISC/README
View File

@ -31,6 +31,7 @@ compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013
compute cnp/atom, Paulo Branicio (USC), branicio at usc.edu, 31 May 2017
compute entropy/atom, Pablo Piaggi (EPFL), pablo.piaggi at epfl.ch, 15 June 2018
compute gyration/shape, Evangelos Voyiatzis, evoyiatzis at gmail.com, 25 July 2019
compute hma, Andrew Schultz & David Kofke (UB), ajs42 at buffalo.edu & kofke at buffalo.edu, 1 Jul 2019
compute pressure/cylinder, Cody K. Addington (NCSU), , 2 Oct 2018
compute momentum, Rupert Nash (University of Edinburgh), r.nash at epcc.ed.ac.uk, 28 June 2019
compute stress/mop, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18

516
src/USER-MISC/compute_hma.cpp Normal file
View File

@ -0,0 +1,516 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
This compute implements harmonically-mapped averaging for crystalline solids.
The current implementation handles atomic crystals.
Computing the heat capacity relies on being able to compute the second
derivative of the energy with respect to atom positions. This capability is
provided by the single2 method in Pair, but is currently only implemented for
the shifted-force LJ potential (lj/smooth/linear). Pair::single2 takes a single
pair and (like Pair::single) returns the energy and sets the force as an out
parameter, but also sets the elements of 6-element double array out parameter,
which are the unique components of the atomic Hessian tensor for the pair. A
helper method exists (Pair::pairTensor), which will compute the tensor from
linear derivatives and the vector between the positions. HMA Heat capacity can
be computed for other models by implementing single2 in those Pair classes.
More information about HMA is available in these publications:
A. J. Schultz, D. A. Kofke, “Comprehensive high-precision high-accuracy
equation of state and coexistence properties for classical Lennard-Jones
crystals and low-temperature fluid phases”, J. Chem. Phys. 149, 204508 (2018)
https://dx.doi.org/10.1063/1.5053714
S. G. Moustafa, A. J. Schultz, D. A. Kofke, “Harmonically Assisted Methods for
Computing the Free Energy of Classical Crystals by Molecular Simulation: A
Comparative Study”, J. Chem. Theory Comput. 13, 825-834 (2017)
https://dx.doi.org/10.1021/acs.jctc.6b01082
S. G. Moustafa, A. J. Schultz, D. A. Kofke, “Very fast averaging of thermal
properties of crystals by molecular simulation”, Phys. Rev. E 92, 043303 (2015)
https://dx.doi.org/10.1103/PhysRevE.92.043303
------------------------------------------------------------------------- */
#include <cmath>
#include <cstring>
#include <mpi.h>
#include "compute_hma.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "group.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "fix_store.h"
#include "memory.h"
#include "error.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "neigh_list.h"
#include <vector>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), id_temp(NULL), deltaR(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal compute hma command");
if (igroup) error->all(FLERR,"Compute hma must use group all");
if (strcmp(arg[3],"NULL") == 0) {error->all(FLERR,"fix ID specifying the set temperature of canonical simulation is required");}
else {
int n = strlen(arg[3]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[3]);
}
create_attribute = 1;
extscalar = 1;
timeflag = 1;
// (from compute displace/atom) create a new fix STORE style
// our new fix's id (id_fix)= compute-ID + COMPUTE_STORE
// our new fix's group = same as compute group
int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
id_fix = new char[n];
strcpy(id_fix,id);
strcat(id_fix,"_COMPUTE_STORE");
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "1";
newarg[5] = (char *) "3";
modify->add_fix(6,newarg);
fix = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;
// calculate xu,yu,zu for fix store array
// skip if reset from restart file
if (fix->restart_reset) fix->restart_reset = 0;
else {
double **xoriginal = fix->astore;
double **x = atom->x;
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
domain->unmap(x[i],image[i],xoriginal[i]);
}
vector_flag = 1;
extvector = -1;
comm_forward = 0; // 3 if 2nd derivative needed
computeU = computeP = computeCv = -1;
returnAnharmonic = 0;
size_vector = 0;
memory->create(extlist, 3, "hma:extlist");
for (int iarg=4; iarg<narg; iarg++) {
if (!strcmp(arg[iarg], "u")) {
if (computeU>-1) continue;
computeU = size_vector;
extlist[size_vector] = 1;
size_vector++;
}
else if (!strcmp(arg[iarg], "p")) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute hma command");
if (computeP>-1) continue;
computeP = size_vector;
deltaPcap = force->numeric(FLERR, arg[iarg+1]);
extlist[size_vector] = 0;
size_vector++;
iarg++;
}
else if (!strcmp(arg[iarg], "cv")) {
if (computeCv>-1) continue;
computeCv = size_vector;
comm_forward = 3;
extlist[size_vector] = 1;
size_vector++;
}
else if (!strcmp(arg[iarg], "anharmonic")) {
// the first time we're called, we'll grab lattice pressure and energy
returnAnharmonic = -1;
}
else {
error->all(FLERR,"Illegal compute hma command");
}
}
if (size_vector == 0) {
error->all(FLERR,"Illegal compute hma command");
}
if (size_vector<3) {
memory->grow(extlist, size_vector, "hma:extlist");
}
memory->create(vector, size_vector, "hma:vector");
if (computeU>-1 || computeCv>-1) {
peflag = 1;
}
if (computeP>-1) {
pressflag = 1;
}
nmax = 0;
}
/* ---------------------------------------------------------------------- */
ComputeHMA::~ComputeHMA()
{
// check nfix in case all fixes have already been deleted
if (modify->nfix) modify->delete_fix(id_fix);
delete [] id_fix;
delete [] id_temp;
memory->destroy(extlist);
memory->destroy(vector);
memory->destroy(deltaR);
}
/* ---------------------------------------------------------------------- */
void ComputeHMA::init() {
if (computeCv>-1) {
if (force->pair == NULL)
error->all(FLERR,"No pair style is defined for compute hma cv");
if (force->pair->single_enable == 0)
error->all(FLERR,"Pair style does not support compute hma cv");
}
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->occasional = 1;
}
void ComputeHMA::init_list(int id, NeighList *ptr)
{
list = ptr;
}
void ComputeHMA::setup()
{
int dummy=0;
int ifix = modify->find_fix(id_temp);
if (ifix < 0) error->all(FLERR,"Could not find compute hma temperature ID");
double * temperat = (double *) modify->fix[ifix]->extract("t_target",dummy);
if (temperat==NULL) error->all(FLERR,"Could not find compute hma temperature ID");
finaltemp = * temperat;
// set fix which stores original atom coords
int ifix2 = modify->find_fix(id_fix);
if (ifix2 < 0) error->all(FLERR,"Could not find hma store fix ID");
fix = (FixStore *) modify->fix[ifix2];
}
/* ---------------------------------------------------------------------- */
void ComputeHMA::compute_vector()
{
invoked_vector = update->ntimestep;
// grow deltaR array if necessary
if (comm_forward>0 && atom->nmax > nmax) {
memory->destroy(deltaR);
nmax = atom->nmax;
memory->create(deltaR,nmax,3,"hma:deltaR");
}
double **xoriginal = fix->astore;
double fdr = 0.0;
double **x = atom->x;
double **f = atom->f;
imageint *image = atom->image;
int nlocal = atom->nlocal;
double *h = domain->h;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double u = 0.0;
if (computeU>-1 || computeCv>-1) {
if (force->pair) u += force->pair->eng_vdwl + force->pair->eng_coul;
if (force->bond) u += force->bond->energy;
if (force->angle) u += force->angle->energy;
if (force->dihedral) u += force->dihedral->energy;
if (force->improper) u += force->improper->energy;
}
int dimension = domain->dimension;
double p = 0, vol = 0;
if (computeP>-1) {
p = virial_compute(dimension);
vol = xprd * yprd;
if (dimension == 3) vol *= zprd;
p *= force->nktv2p / (dimension*vol);
if (returnAnharmonic == -1) {
pLat = p;
}
}
if (domain->triclinic == 0) {
for (int i = 0; i < nlocal; i++) {
int xbox = (image[i] & IMGMASK) - IMGMAX;
int ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
int zbox = (image[i] >> IMG2BITS) - IMGMAX;
double dx = x[i][0] + xbox*xprd - xoriginal[i][0];
double dy = x[i][1] + ybox*yprd - xoriginal[i][1];
double dz = x[i][2] + zbox*zprd - xoriginal[i][2];
if (comm_forward>0) {
deltaR[i][0] = dx;
deltaR[i][1] = dy;
deltaR[i][2] = dz;
}
fdr += dx*f[i][0] + dy*f[i][1] + dz*f[i][2];
}
}
else {
for (int i = 0; i < nlocal; i++) {
int xbox = (image[i] & IMGMASK) - IMGMAX;
int ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
int zbox = (image[i] >> IMG2BITS) - IMGMAX;
double dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox - xoriginal[i][0];
double dy = x[i][1] + h[1]*ybox + h[3]*zbox - xoriginal[i][1];
double dz = x[i][2] + h[2]*zbox - xoriginal[i][2];
if (comm_forward>0) {
deltaR[i][0] = dx;
deltaR[i][1] = dy;
deltaR[i][2] = dz;
}
fdr += ((dx*f[i][0])+(dy*f[i][1])+(dz*f[i][2]));
}
}
double phiSum = 0.0;
if (computeCv>-1) {
comm->forward_comm_compute(this);
int *type = atom->type;
double** cutsq = force->pair->cutsq;
if (force->pair) {
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
double *special_coul = force->special_coul;
int newton_pair = force->newton_pair;
if (update->firststep == update->ntimestep) neighbor->build_one(list,1);
else neighbor->build_one(list);
int inum = list->inum;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
for (int ii = 0; ii < inum; ii++) {
int i = ilist[ii];
double fac = (newton_pair || i < nlocal) ? 1.0 : 0.5;
double* ix = x[i];
int itype = type[i];
int *jlist = firstneigh[i];
int jnum = numneigh[i];
double *idr = deltaR[i];
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
if (!newton_pair && j>=nlocal) fac -= 0.5;
double factor_lj = special_lj[sbmask(j)];
double factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
double* jx = x[j];
double delr[3];
delr[0] = ix[0] - jx[0];
delr[1] = ix[1] - jx[1];
delr[2] = ix[2] - jx[2];
double rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
int jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
double* jdr = deltaR[j];
double fforce, d2u[6];
force->pair->single_hessian(i, j, itype, jtype, rsq, delr, factor_coul, factor_lj, fforce, d2u);
int m = 0;
for (int k=0; k<3; k++) {
double a = fac;
for (int l=k; l<3; l++) {
phiSum += a*(idr[k]*jdr[l]+jdr[k]*idr[l])*d2u[m];
phiSum -= a*(idr[k]*idr[l]*d2u[m] + jdr[k]*jdr[l]*d2u[m]);
m++;
if (k==l) a *= 2;
}
}
}
}
}
}
}
// compute and sum up properties across processors
double fdrTotal;
MPI_Allreduce(&fdr,&fdrTotal,1,MPI_DOUBLE,MPI_SUM,world);
double uTotal;
if (computeU>-1 || computeCv>-1) {
MPI_Allreduce(&u,&uTotal,1,MPI_DOUBLE,MPI_SUM,world);
if (returnAnharmonic == -1) {
uLat = uTotal;
}
if (computeU>-1) {
if (returnAnharmonic) {
vector[computeU] = uTotal - uLat + 0.5*fdrTotal;
}
else {
vector[computeU] = uTotal + 0.5*fdrTotal + 0.5*dimension*(atom->natoms - 1)*force->boltz*finaltemp;
}
}
}
if (computeP>-1) {
double fv = ((deltaPcap)-(force->boltz*finaltemp*force->nktv2p*atom->natoms/vol))/(force->boltz*finaltemp*dimension*(atom->natoms - 1));
if (returnAnharmonic) {
vector[computeP] = p - pLat + (fv*fdrTotal);
}
else {
vector[computeP] = p + (fv*fdrTotal) + deltaPcap;
}
}
if (computeCv>-1) {
if (computeU==-1) MPI_Allreduce(&u,&uTotal,1,MPI_DOUBLE,MPI_SUM,world);
double buTot;
if (returnAnharmonic) {
buTot = (uTotal - uLat + 0.5*fdrTotal)/finaltemp;
}
else {
buTot = (uTotal + 0.5*fdrTotal)/finaltemp + 0.5*dimension*(atom->natoms - 1)*force->boltz;
}
double one = -0.25*(fdr + phiSum)/finaltemp;
double Cv;
MPI_Allreduce(&one,&Cv,1,MPI_DOUBLE,MPI_SUM,world);
vector[computeCv] = Cv + buTot*buTot;
if (!returnAnharmonic) {
vector[computeCv] += 0.5*dimension*(atom->natoms-1);
}
}
if (returnAnharmonic == -1) {
// we have our lattice properties
returnAnharmonic = 1;
}
}
double ComputeHMA::virial_compute(int n)
{
double v = 0;
// sum contributions to virial from forces and fixes
if (force->pair) v += sumVirial(n, force->pair->virial);
if (force->bond) v += sumVirial(n, force->bond->virial);
if (force->angle) v += sumVirial(n, force->angle->virial);
if (force->dihedral) v += sumVirial(n, force->dihedral->virial);
if (force->improper) v += sumVirial(n, force->improper->virial);
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->thermo_virial) v += sumVirial(n, modify->fix[i]->virial);
// sum virial across procs
double virial;
MPI_Allreduce(&v,&virial,1,MPI_DOUBLE,MPI_SUM,world);
// KSpace virial contribution is already summed across procs
if (force->kspace)
for (int i = 0; i < n; i++) virial += force->kspace->virial[i];
return virial;
}
/* ---------------------------------------------------------------------- */
int ComputeHMA::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
double **xoriginal = fix->astore;
imageint *image = atom->image;
double **x = atom->x;
double *h = domain->h;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
int m = 0;
for (int ii = 0; ii < n; ii++) {
int i = list[ii];
buf[m++] = deltaR[i][0];
buf[m++] = deltaR[i][1];
buf[m++] = deltaR[i][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeHMA::unpack_forward_comm(int n, int first, double *buf)
{
double **xoriginal = fix->astore;
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
deltaR[i][0] = buf[m++];
deltaR[i][1] = buf[m++];
deltaR[i][2] = buf[m++];
}
}
/* ----------------------------------------------------------------------
initialize one atom's storage values, called when atom is created
------------------------------------------------------------------------- */
void ComputeHMA::set_arrays(int i)
{
double **xoriginal = fix->astore;
double **x = atom->x;
xoriginal[i][0] = x[i][0];
xoriginal[i][1] = x[i][1];
xoriginal[i][2] = x[i][2];
}
double ComputeHMA::memory_usage()
{
double bytes = nmax * 3 * sizeof(double);
return bytes;
}

67
src/USER-MISC/compute_hma.h Normal file
View File

@ -0,0 +1,67 @@
/* -*- 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 COMPUTE_CLASS
ComputeStyle(hma,ComputeHMA)
#else
#ifndef LMP_COMPUTE_HMA_H
#define LMP_COMPUTE_HMA_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeHMA : public Compute {
public:
ComputeHMA(class LAMMPS *, int, char **);
~ComputeHMA();
void setup();
void init();
void init_list(int, class NeighList *);
void compute_vector();
void set_arrays(int);
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
double memory_usage();
private:
int nmax;
int atomsingroup;
char *id_fix;
char *id_temp;
double finaltemp;
class FixStore *fix;
double boltz, nktv2p, inv_volume;
double deltaPcap;
double virial_compute(int);
static double sumVirial(int n, double* v) {
double x = 0;
for (int i=0; i<n; i++) x += v[i];
return x;
}
int computeU, computeP, computeCv;
class NeighList *list; // half neighbor list
double **deltaR;
int returnAnharmonic;
double uLat, pLat;
};
}
#endif
#endif

13
src/pair.cpp
View File

@ -54,6 +54,7 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp)
comm_forward = comm_reverse = comm_reverse_off = 0;
single_enable = 1;
single_hessian_enable = 0;
restartinfo = 1;
respa_enable = 0;
one_coeff = 0;
@ -1743,6 +1744,18 @@ void Pair::init_bitmap(double inner, double outer, int ntablebits,
/* ---------------------------------------------------------------------- */
void Pair::hessian_twobody(double fforce, double dfac, double delr[3], double phiTensor[6]) {
int m = 0;
for (int k=0; k<3; k++) {
phiTensor[m] = fforce;
for (int l=k; l<3; l++) {
if (l>k) phiTensor[m] = 0;
phiTensor[m++] += delr[k]*delr[l] * dfac;
}
}
}
/* ---------------------------------------------------------------------- */
double Pair::memory_usage()
{
double bytes = comm->nthreads*maxeatom * sizeof(double);

11
src/pair.h
View File

@ -46,6 +46,7 @@ class Pair : protected Pointers {
int comm_reverse_off; // size of reverse comm even if newton off
int single_enable; // 1 if single() routine exists
int single_hessian_enable; // 1 if single_hessian() routine exists
int restartinfo; // 1 if pair style writes restart info
int respa_enable; // 1 if inner/middle/outer rRESPA routines
int one_coeff; // 1 if allows only one coeff * * call
@ -148,6 +149,16 @@ class Pair : protected Pointers {
return 0.0;
}
void hessian_twobody(double fforce, double dfac, double delr[3], double phiTensor[6]);
virtual double single_hessian(int, int, int, int,
double, double[3], double, double,
double& fforce, double d2u[6]) {
fforce = 0.0;
for (int i=0; i<6; i++) d2u[i] = 0;
return 0.0;
}
virtual void settings(int, char **) = 0;
virtual void coeff(int, char **) = 0;

27
src/pair_lj_smooth_linear.cpp
View File

@ -29,7 +29,9 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairLJSmoothLinear::PairLJSmoothLinear(LAMMPS *lmp) : Pair(lmp) {}
PairLJSmoothLinear::PairLJSmoothLinear(LAMMPS *lmp) : Pair(lmp) {
single_hessian_enable = 1;
}
/* ---------------------------------------------------------------------- */
@ -345,3 +347,26 @@ double PairLJSmoothLinear::single(int /*i*/, int /*j*/, int itype, int jtype,
return factor_lj*philj;
}
double PairLJSmoothLinear::single_hessian(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
double delr[3], double /*factor_coul*/, double factor_lj,
double &fforce, double d2u[6])
{
double r2inv,r6inv,forcelj,philj,r,rinv;
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
rinv = sqrt(r2inv);
r = sqrt(rsq);
forcelj = r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
forcelj = rinv*forcelj - dljcut[itype][jtype];
fforce = factor_lj*forcelj*rinv;
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
philj = philj - ljcut[itype][jtype]
+ (r-cut[itype][jtype])*dljcut[itype][jtype];
double d2r = factor_lj * r6inv * (13.0*lj1[itype][jtype]*r6inv - 7.0*lj2[itype][jtype])/rsq;
hessian_twobody(fforce, -(fforce + d2r) / rsq, delr, d2u);
return factor_lj*philj;
}

1
src/pair_lj_smooth_linear.h
View File

@ -38,6 +38,7 @@ class PairLJSmoothLinear : public Pair {
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
double single_hessian(int, int, int, int, double, double[3], double, double, double&, double[6]);
protected:
double cut_global;