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

Merge remote-tracking branch 'origin/master' into doc-reorg-intro-howto

This commit is contained in:
Richard Berger 2018-08-03 00:12:58 -04:00
parent 9cc75792d6 1b0a8fdc9b
commit 2321c8ff37
27 changed files with 898 additions and 782 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.github/CODEOWNERS vendored
View File

@ -45,6 +45,7 @@ src/USER-MISC/*_grem.* @dstelter92
# tools
tools/msi2lmp/* @akohlmey
tools/emacs/* @HaoZeke
# cmake
cmake/* @junghans @rbberger

1
cmake/CMakeLists.txt
View File

@ -181,7 +181,6 @@ endmacro()
pkg_depends(MPIIO MPI)
pkg_depends(USER-ATC MANYBODY)
pkg_depends(USER-LB MPI)
pkg_depends(USER-MISC MANYBODY)
pkg_depends(USER-PHONON KSPACE)
######################################################

4
doc/src/Howto_body.txt
View File

@ -151,8 +151,8 @@ center-of-mass position of the particle is specified by the x,y,z
values in the {Atoms} section of the data file, as is the total mass
of the body particle.
The "pair_style body"_pair_body.html command can be used with this
body style to compute body/body and body/non-body interactions.
The "pair_style body/nparticle"_pair_body_nparticle.html command can be used
with this body style to compute body/body and body/non-body interactions.
For output purposes via the "compute
body/local"_compute_body_local.html and "dump local"_dump.html

4
doc/src/Howto_spherical.txt
View File

@ -129,7 +129,7 @@ styles"_pair_style.html that generate torque:
"pair_style lubricate"_pair_lubricate.html
"pair_style line/lj"_pair_line_lj.html
"pair_style tri/lj"_pair_tri_lj.html
"pair_style body"_pair_body.html :ul
"pair_style body/nparticle"_pair_body_nparticle.html :ul
The granular pair styles are used with spherical particles. The
dipole pair style is used with the dipole atom style, which could be
@ -240,4 +240,4 @@ list of sub-particles. Individual body partices are typically treated
as rigid bodies, and their motion integrated with a command like "fix
nve/body"_fix_nve_body.html. Interactions between pairs of body
particles are computed via a command like "pair_style
body"_pair_body.html.
body/nparticle"_pair_body_nparticle.html.

4
doc/src/Manual.txt
View File

@ -1,7 +1,7 @@
<!-- HTML_ONLY -->
<HEAD>
<TITLE>LAMMPS Users Manual</TITLE>
<META NAME="docnumber" CONTENT="16 Jul 2018 version">
<META NAME="docnumber" CONTENT="2 Aug 2018 version">
<META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
<META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation. This software and manual is distributed under the GNU General Public License.">
</HEAD>
@ -21,7 +21,7 @@
<H1></H1>
LAMMPS Documentation :c,h1
16 Jul 2018 version :c,h2
2 Aug 2018 version :c,h2
"What is a LAMMPS version?"_Manual_version.html

2
doc/src/Packages_details.txt
View File

@ -148,7 +148,7 @@ src/BODY filenames -> commands
"Howto_body"_Howto_body.html
"atom_style body"_atom_style.html
"fix nve/body"_fix_nve_body.html
"pair_style body"_pair_body.html
"pair_style body/nparticle"_pair_body_nparticle.html
examples/body :ul
:line

2
doc/src/Section_commands.txt
View File

@ -969,6 +969,7 @@ OPT.
"dsmc"_pair_dsmc.html,
"eam (gikot)"_pair_eam.html,
"eam/alloy (gikot)"_pair_eam.html,
"eam/cd (o)"_pair_eam.html,
"eam/fs (gikot)"_pair_eam.html,
"eim (o)"_pair_eim.html,
"gauss (go)"_pair_gauss.html,
@ -1069,7 +1070,6 @@ package"_Section_start.html#start_3.
"coul/shield"_pair_coul_shield.html,
"dpd/fdt"_pair_dpd_fdt.html,
"dpd/fdt/energy (k)"_pair_dpd_fdt.html,
"eam/cd (o)"_pair_eam.html,
"edip (o)"_pair_edip.html,
"edip/multi"_pair_edip.html,
"edpd"_pair_meso.html,

6
doc/src/Tools.txt
View File

@ -249,9 +249,9 @@ These tools were provided by Andres Jaramillo-Botero at CalTech
emacs tool :h4,link(emacs)
The tools/emacs directory contains a Lips add-on file for Emacs that
enables a lammps-mode for editing of input scripts when using Emacs,
with various highlighting options setup.
The tools/emacs directory contains an Emacs Lisp add-on file for GNU Emacs
that enables a lammps-mode for editing input scripts when using GNU Emacs,
with various highlighting options set up.
These tools were provided by Aidan Thompson at Sandia
(athomps at sandia.gov).

2
doc/src/body.txt
View File

@ -151,7 +151,7 @@ center-of-mass position of the particle is specified by the x,y,z
values in the {Atoms} section of the data file, as is the total mass
of the body particle.
The "pair_style body"_pair_body.html command can be used with this
The "pair_style body/nparticle"_pair_body_nparticle.html command can be used with this
body style to compute body/body and body/non-body interactions.
For output purposes via the "compute

2
doc/src/fixes.txt
View File

@ -168,6 +168,8 @@ Fixes :h1
fix_viscosity
fix_viscous
fix_wall
fix_wall_body_polygon
fix_wall_body_polyhedron
fix_wall_ees
fix_wall_gran
fix_wall_gran_region

9
doc/src/pair_eam.txt
View File

@ -413,15 +413,10 @@ The eam pair styles can only be used via the {pair} keyword of the
[Restrictions:]
All of these styles except the {eam/cd} style are part of the MANYBODY
package. They are only enabled if LAMMPS was built with that package.
All of these styles are part of the MANYBODY package. They are only
enabled if LAMMPS was built with that package.
See the "Making LAMMPS"_Section_start.html#start_3 section for more info.
The {eam/cd} style is part of the USER-MISC package and also requires
the MANYBODY package. It is only enabled if LAMMPS was built with
those packages. See the "Making LAMMPS"_Section_start.html#start_3
section for more info.
[Related commands:]
"pair_coeff"_pair_coeff.html

2
doc/src/pair_style.txt
View File

@ -104,7 +104,7 @@ in the pair section of "this page"_Section_commands.html#cmd_5.
"pair_style airebo"_pair_airebo.html - AIREBO potential of Stuart
"pair_style airebo/morse"_pair_airebo.html - AIREBO with Morse instead of LJ
"pair_style beck"_pair_beck.html - Beck potential
"pair_style body"_pair_body.html - interactions between body particles
"pair_style body/nparticle"_pair_body_nparticle.html - interactions between body particles
"pair_style bop"_pair_bop.html - BOP potential of Pettifor
"pair_style born"_pair_born.html - Born-Mayer-Huggins potential
"pair_style born/coul/long"_pair_born.html - Born-Mayer-Huggins with long-range Coulombics

3
doc/src/pairs.txt
View File

@ -10,8 +10,9 @@ Pair Styles :h1
pair_airebo
pair_awpmd
pair_beck
pair_body
pair_body_nparticle
pair_body_rounded_polygon
pair_body_rounded_polyhedron
pair_bop
pair_born
pair_brownian

2
examples/USER/diffraction/BulkNi.in
View File

@ -20,7 +20,7 @@ compute XRD all xrd 1.541838 Ni 2Theta 40 80 c 2 2 2 LP 1 echo
compute SAED all saed 0.0251 Ni Kmax 0.85 Zone 1 0 0 c 0.025 0.025 0.025 &
dR_Ewald 0.05 echo manual
fix 1 all ave/histo 1 1 1 40 80 200 c_XRD[1] c_XRD[2] &
fix 1 all ave/histo/weight 1 1 1 40 80 200 c_XRD[1] c_XRD[2] &
mode vector file $A.hist.xrd
fix 2 all saed/vtk 1 1 1 c_SAED file $A_001.saed

2
src/KOKKOS/fix_nh_kokkos.cpp
View File

@ -148,7 +148,7 @@ void FixNHKokkos<DeviceType>::setup(int vflag)
if (pstat_flag) {
double kt = boltz * t_target;
double nkt = atom->natoms * kt;
double nkt = (atom->natoms + 1) * kt;
for (int i = 0; i < 3; i++)
if (p_flag[i])

677
src/MANYBODY/pair_eam_cd.cpp Normal file
View File

@ -0,0 +1,677 @@
/* ----------------------------------------------------------------------
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: Alexander Stukowski
Technical University of Darmstadt,
Germany Department of Materials Science
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "pair_eam_cd.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define ASSERT(cond)
#define MAXLINE 1024 // This sets the maximum line length in EAM input files.
PairEAMCD::PairEAMCD(LAMMPS *lmp, int _cdeamVersion)
: PairEAM(lmp), PairEAMAlloy(lmp), cdeamVersion(_cdeamVersion)
{
single_enable = 0;
restartinfo = 0;
rhoB = NULL;
D_values = NULL;
hcoeff = NULL;
// Set communication buffer sizes needed by this pair style.
if (cdeamVersion == 1) {
comm_forward = 4;
comm_reverse = 3;
} else if (cdeamVersion == 2) {
comm_forward = 3;
comm_reverse = 2;
} else {
error->all(FLERR,"Invalid eam/cd potential version.");
}
}
PairEAMCD::~PairEAMCD()
{
memory->destroy(rhoB);
memory->destroy(D_values);
if (hcoeff) delete[] hcoeff;
}
void PairEAMCD::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,rhoip,rhojp,recip,phi;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
// Grow per-atom arrays if necessary
if (atom->nmax > nmax) {
memory->destroy(rho);
memory->destroy(fp);
memory->destroy(rhoB);
memory->destroy(D_values);
nmax = atom->nmax;
memory->create(rho,nmax,"pair:rho");
memory->create(rhoB,nmax,"pair:rhoB");
memory->create(fp,nmax,"pair:fp");
memory->create(D_values,nmax,"pair:D_values");
}
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// Zero out per-atom arrays.
int m = nlocal + atom->nghost;
for (i = 0; i < m; i++) {
rho[i] = 0.0;
rhoB[i] = 0.0;
D_values[i] = 0.0;
}
// Stage I
// Compute rho and rhoB at each local atom site.
// Additionally calculate the D_i values here if we are using the
// one-site formulation. For the two-site formulation we have to
// calculate the D values in an extra loop (Stage II).
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
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];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutforcesq) {
jtype = type[j];
double r = sqrt(rsq);
const EAMTableIndex index = radiusToTableIndex(r);
double localrho = RhoOfR(index, jtype, itype);
rho[i] += localrho;
if (jtype == speciesB) rhoB[i] += localrho;
if (newton_pair || j < nlocal) {
localrho = RhoOfR(index, itype, jtype);
rho[j] += localrho;
if (itype == speciesB) rhoB[j] += localrho;
}
if (cdeamVersion == 1 && itype != jtype) {
// Note: if the i-j interaction is not concentration dependent (because either
// i or j are not species A or B) then its contribution to D_i and D_j should
// be ignored.
// This if-clause is only required for a ternary.
if ((itype == speciesA && jtype == speciesB)
|| (jtype == speciesA && itype == speciesB)) {
double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);
D_values[i] += Phi_AB;
if (newton_pair || j < nlocal)
D_values[j] += Phi_AB;
}
}
}
}
}
// Communicate and sum densities.
if (newton_pair) {
communicationStage = 1;
comm->reverse_comm_pair(this);
}
// fp = derivative of embedding energy at each atom
// phi = embedding energy at each atom
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
EAMTableIndex index = rhoToTableIndex(rho[i]);
fp[i] = FPrimeOfRho(index, type[i]);
if (eflag) {
phi = FofRho(index, type[i]);
if (eflag_global) eng_vdwl += phi;
if (eflag_atom) eatom[i] += phi;
}
}
// Communicate derivative of embedding function and densities
// and D_values (this for one-site formulation only).
communicationStage = 2;
comm->forward_comm_pair(this);
// The electron densities may not drop to zero because then the
// concentration would no longer be defined. But the concentration
// is not needed anyway if there is no interaction with another atom,
// which is the case if the electron density is exactly zero.
// That's why the following lines have been commented out.
//
//for (i = 0; i < nlocal + atom->nghost; i++) {
// if (rho[i] == 0 && (type[i] == speciesA || type[i] == speciesB))
// error->one(FLERR,"CD-EAM potential routine: Detected atom with zero electron density.");
//}
// Stage II
// This is only required for the original two-site formulation of the CD-EAM potential.
if (cdeamVersion == 2) {
// Compute intermediate value D_i for each atom.
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
// This code line is required for ternary alloys.
if (itype != speciesA && itype != speciesB) continue;
double x_i = rhoB[i] / rho[i]; // Concentration at atom i.
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
if (itype == jtype) continue;
// This code line is required for ternary alloys.
if (jtype != speciesA && jtype != speciesB) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutforcesq) {
double r = sqrt(rsq);
const EAMTableIndex index = radiusToTableIndex(r);
// The concentration independent part of the cross pair potential.
double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);
// Average concentration of two sites
double x_ij = 0.5 * (x_i + rhoB[j]/rho[j]);
// Calculate derivative of h(x_ij) polynomial function.
double h_prime = evalHprime(x_ij);
D_values[i] += h_prime * Phi_AB / (2.0 * rho[i] * rho[i]);
if (newton_pair || j < nlocal)
D_values[j] += h_prime * Phi_AB / (2.0 * rho[j] * rho[j]);
}
}
}
// Communicate and sum D values.
if (newton_pair) {
communicationStage = 3;
comm->reverse_comm_pair(this);
}
communicationStage = 4;
comm->forward_comm_pair(this);
}
// Stage III
// Compute force acting on each atom.
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
// Concentration at site i
// The value -1 indicates: no concentration dependence for all interactions of atom i.
// It will be replaced by the concentration at site i if atom i is either A or B.
double x_i = -1.0;
double D_i, h_prime_i;
// This if-clause is only required for ternary alloys.
if ((itype == speciesA || itype == speciesB) && rho[i] != 0.0) {
// Compute local concentration at site i.
x_i = rhoB[i]/rho[i];
ASSERT(x_i >= 0 && x_i<=1.0);
if (cdeamVersion == 1) {
// Calculate derivative of h(x_i) polynomial function.
h_prime_i = evalHprime(x_i);
D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]);
} else if (cdeamVersion == 2) {
D_i = D_values[i];
} else {
ASSERT(false);
}
}
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutforcesq) {
jtype = type[j];
double r = sqrt(rsq);
const EAMTableIndex index = radiusToTableIndex(r);
// rhoip = derivative of (density at atom j due to atom i)
// rhojp = derivative of (density at atom i due to atom j)
// psip needs both fp[i] and fp[j] terms since r_ij appears in two
// terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)
// hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
rhoip = RhoPrimeOfR(index, itype, jtype);
rhojp = RhoPrimeOfR(index, jtype, itype);
fpair = fp[i]*rhojp + fp[j]*rhoip;
recip = 1.0/r;
// The value -1 indicates: no concentration dependence for this
// i-j pair because atom j is not of species A nor B.
double x_j = -1;
// This code line is required for ternary alloy.
if (jtype == speciesA || jtype == speciesB) {
ASSERT(rho[i] != 0.0);
ASSERT(rho[j] != 0.0);
// Compute local concentration at site j.
x_j = rhoB[j]/rho[j];
ASSERT(x_j >= 0 && x_j<=1.0);
double D_j=0.0;
if (cdeamVersion == 1) {
// Calculate derivative of h(x_j) polynomial function.
double h_prime_j = evalHprime(x_j);
D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]);
} else if (cdeamVersion == 2) {
D_j = D_values[j];
} else {
ASSERT(false);
}
double t2 = -rhoB[j];
if (itype == speciesB) t2 += rho[j];
fpair += D_j * rhoip * t2;
}
// This if-clause is only required for a ternary alloy.
// Actually we don't need it at all because D_i should be zero
// anyway if atom i has no concentration dependent interactions
// (because it is not species A or B).
if (x_i != -1.0) {
double t1 = -rhoB[i];
if (jtype == speciesB) t1 += rho[i];
fpair += D_i * rhojp * t1;
}
double phip;
double phi = PhiOfR(index, itype, jtype, recip, phip);
if (itype == jtype || x_i == -1.0 || x_j == -1.0) {
// Case of no concentration dependence.
fpair += phip;
} else {
// We have a concentration dependence for the i-j interaction.
double h=0.0;
if (cdeamVersion == 1) {
// Calculate h(x_i) polynomial function.
double h_i = evalH(x_i);
// Calculate h(x_j) polynomial function.
double h_j = evalH(x_j);
h = 0.5 * (h_i + h_j);
} else if (cdeamVersion == 2) {
// Average concentration.
double x_ij = 0.5 * (x_i + x_j);
// Calculate h(x_ij) polynomial function.
h = evalH(x_ij);
} else {
ASSERT(false);
}
fpair += h * phip;
phi *= h;
}
// Divide by r_ij and negate to get forces from gradient.
fpair /= -r;
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) evdwl = phi;
if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairEAMCD::coeff(int narg, char **arg)
{
PairEAMAlloy::coeff(narg, arg);
// Make sure the EAM file is a CD-EAM binary alloy.
if (setfl->nelements < 2)
error->all(FLERR,"The EAM file must contain at least 2 elements to be used with the eam/cd pair style.");
// Read in the coefficients of the h polynomial from the end of the EAM file.
read_h_coeff(arg[2]);
// Determine which atom type is the A species and which is the B
// species in the alloy. By default take the first element (index 0)
// in the EAM file as the A species and the second element (index 1)
// in the EAM file as the B species.
speciesA = -1;
speciesB = -1;
for (int i = 1; i <= atom->ntypes; i++) {
if (map[i] == 0) {
if (speciesA >= 0)
error->all(FLERR,"The first element from the EAM file may only be mapped to a single atom type.");
speciesA = i;
}
if (map[i] == 1) {
if (speciesB >= 0)
error->all(FLERR,"The second element from the EAM file may only be mapped to a single atom type.");
speciesB = i;
}
}
if (speciesA < 0)
error->all(FLERR,"The first element from the EAM file must be mapped to exactly one atom type.");
if (speciesB < 0)
error->all(FLERR,"The second element from the EAM file must be mapped to exactly one atom type.");
}
/* ----------------------------------------------------------------------
Reads in the h(x) polynomial coefficients
------------------------------------------------------------------------- */
void PairEAMCD::read_h_coeff(char *filename)
{
if (comm->me == 0) {
// Open potential file
FILE *fptr;
char line[MAXLINE];
char nextline[MAXLINE];
fptr = force->open_potential(filename);
if (fptr == NULL) {
char str[128];
sprintf(str,"Cannot open EAM potential file %s", filename);
error->one(FLERR,str);
}
// h coefficients are stored at the end of the file.
// Skip to last line of file.
while(fgets(nextline, MAXLINE, fptr) != NULL) {
strcpy(line, nextline);
}
char* ptr = strtok(line, " \t\n\r\f");
int degree = atoi(ptr);
nhcoeff = degree+1;
hcoeff = new double[nhcoeff];
int i = 0;
while((ptr = strtok(NULL," \t\n\r\f")) != NULL && i < nhcoeff) {
hcoeff[i++] = atof(ptr);
}
if (i != nhcoeff || nhcoeff < 1)
error->one(FLERR,"Failed to read h(x) function coefficients from EAM file.");
// Close the potential file.
fclose(fptr);
}
MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world);
if (comm->me != 0) hcoeff = new double[nhcoeff];
MPI_Bcast(hcoeff, nhcoeff, MPI_DOUBLE, 0, world);
}
/* ---------------------------------------------------------------------- */
int PairEAMCD::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
m = 0;
if (communicationStage == 2) {
if (cdeamVersion == 1) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = fp[j];
buf[m++] = rho[j];
buf[m++] = rhoB[j];
buf[m++] = D_values[j];
}
return m;
} else if (cdeamVersion == 2) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = fp[j];
buf[m++] = rho[j];
buf[m++] = rhoB[j];
}
return m;
} else { ASSERT(false); return 0; }
} else if (communicationStage == 4) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = D_values[j];
}
return m;
} else return 0;
}
/* ---------------------------------------------------------------------- */
void PairEAMCD::unpack_forward_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if (communicationStage == 2) {
if (cdeamVersion == 1) {
for (i = first; i < last; i++) {
fp[i] = buf[m++];
rho[i] = buf[m++];
rhoB[i] = buf[m++];
D_values[i] = buf[m++];
}
} else if (cdeamVersion == 2) {
for (i = first; i < last; i++) {
fp[i] = buf[m++];
rho[i] = buf[m++];
rhoB[i] = buf[m++];
}
} else {
ASSERT(false);
}
} else if (communicationStage == 4) {
for (i = first; i < last; i++) {
D_values[i] = buf[m++];
}
}
}
/* ---------------------------------------------------------------------- */
int PairEAMCD::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if (communicationStage == 1) {
if (cdeamVersion == 1) {
for (i = first; i < last; i++) {
buf[m++] = rho[i];
buf[m++] = rhoB[i];
buf[m++] = D_values[i];
}
return m;
} else if (cdeamVersion == 2) {
for (i = first; i < last; i++) {
buf[m++] = rho[i];
buf[m++] = rhoB[i];
}
return m;
} else { ASSERT(false); return 0; }
} else if (communicationStage == 3) {
for (i = first; i < last; i++) {
buf[m++] = D_values[i];
}
return m;
} else return 0;
}
/* ---------------------------------------------------------------------- */
void PairEAMCD::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
if (communicationStage == 1) {
if (cdeamVersion == 1) {
for (i = 0; i < n; i++) {
j = list[i];
rho[j] += buf[m++];
rhoB[j] += buf[m++];
D_values[j] += buf[m++];
}
} else if (cdeamVersion == 2) {
for (i = 0; i < n; i++) {
j = list[i];
rho[j] += buf[m++];
rhoB[j] += buf[m++];
}
} else {
ASSERT(false);
}
} else if (communicationStage == 3) {
for (i = 0; i < n; i++) {
j = list[i];
D_values[j] += buf[m++];
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairEAMCD::memory_usage()
{
double bytes = 2 * nmax * sizeof(double);
return PairEAMAlloy::memory_usage() + bytes;
}

22
src/USER-MISC/pair_cdeam.h → src/MANYBODY/pair_eam_cd.h
View File

@ -13,26 +13,26 @@
#ifdef PAIR_CLASS
PairStyle(eam/cd,PairCDEAM_OneSite)
PairStyle(eam/cd/old,PairCDEAM_TwoSite)
PairStyle(eam/cd,PairEAMCD_OneSite)
PairStyle(eam/cd/old,PairEAMCD_TwoSite)
#else
#ifndef LMP_PAIR_CDEAM_H
#define LMP_PAIR_CDEAM_H
#ifndef LMP_PAIR_EAM_CD_H
#define LMP_PAIR_EAM_CD_H
#include "pair_eam_alloy.h"
namespace LAMMPS_NS {
class PairCDEAM : public PairEAMAlloy
class PairEAMCD : public PairEAMAlloy
{
public:
/// Constructor.
PairCDEAM(class LAMMPS*, int cdeamVersion);
PairEAMCD(class LAMMPS*, int cdeamVersion);
/// Destructor.
virtual ~PairCDEAM();
virtual ~PairEAMCD();
/// Calculates the energies and forces for all atoms in the system.
virtual void compute(int, int);
@ -211,19 +211,19 @@ public:
};
/// The one-site concentration formulation of CD-EAM.
class PairCDEAM_OneSite : public PairCDEAM
class PairEAMCD_OneSite : public PairEAMCD
{
public:
/// Constructor.
PairCDEAM_OneSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 1) {}
PairEAMCD_OneSite(class LAMMPS* lmp) : PairEAM(lmp), PairEAMCD(lmp, 1) {}
};
/// The two-site concentration formulation of CD-EAM.
class PairCDEAM_TwoSite : public PairCDEAM
class PairEAMCD_TwoSite : public PairEAMCD
{
public:
/// Constructor.
PairCDEAM_TwoSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 2) {}
PairEAMCD_TwoSite(class LAMMPS* lmp) : PairEAM(lmp), PairEAMCD(lmp, 2) {}
};
}

3
src/Purge.list
View File

@ -24,6 +24,9 @@ style_nstencil.h
style_ntopo.h
# other auto-generated files
lmpinstalledpkgs.h
# renamed on 31 July 2018
pair_cdeam.h
pair_cdeam.cpp
# renamed on 20 July 2018
pair_body.h
pair_body.cpp

45
src/USER-BOCS/fix_bocs.cpp
View File

@ -846,7 +846,7 @@ void FixBocs::setup(int vflag)
if (pstat_flag) {
double kt = boltz * t_target;
double nkt = atom->natoms * kt;
double nkt = (atom->natoms + 1) * kt;
for (int i = 0; i < 3; i++)
if (p_flag[i])
@ -1508,7 +1508,7 @@ double FixBocs::compute_scalar()
double volume;
double energy;
double kt = boltz * t_target;
double lkt_press = kt;
double lkt_press = 0.0;
int ich;
if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
else volume = domain->xprd * domain->yprd;
@ -1539,15 +1539,21 @@ double FixBocs::compute_scalar()
// sum is over barostatted dimensions
if (pstat_flag) {
for (i = 0; i < 3; i++)
if (p_flag[i])
for (i = 0; i < 3; i++) {
if (p_flag[i]) {
energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i] +
p_hydro*(volume-vol0) / (pdim*nktv2p);
lkt_press += kt;
}
}
if (pstyle == TRICLINIC) {
for (i = 3; i < 6; i++)
if (p_flag[i])
for (i = 3; i < 6; i++) {
if (p_flag[i]) {
energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i];
lkt_press += kt;
}
}
}
// extra contributions from thermostat chain for barostat
@ -1880,15 +1886,14 @@ void FixBocs::nhc_temp_integrate()
void FixBocs::nhc_press_integrate()
{
int ich,i;
int ich,i,pdof;
double expfac,factor_etap,kecurrent;
double kt = boltz * t_target;
double lkt_press = kt;
// Update masses, to preserve initial freq, if flag set
if (omega_mass_flag) {
double nkt = atom->natoms * kt;
double nkt = (atom->natoms + 1) * kt;
for (int i = 0; i < 3; i++)
if (p_flag[i])
omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
@ -1912,14 +1917,24 @@ void FixBocs::nhc_press_integrate()
}
kecurrent = 0.0;
for (i = 0; i < 3; i++)
if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
if (pstyle == TRICLINIC) {
for (i = 3; i < 6; i++)
if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
pdof = 0;
for (i = 0; i < 3; i++) {
if (p_flag[i]) {
kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
pdof++;
}
}
if (pstyle == TRICLINIC) {
for (i = 3; i < 6; i++) {
if (p_flag[i]) {
kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
pdof++;
}
}
}
double lkt_press = pdof * kt;
etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
double ncfac = 1.0/nc_pchain;

40
src/USER-MISC/Install.sh
View File

@ -1,40 +0,0 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
mode=$1
# enforce using portable C locale
LC_ALL=C
export LC_ALL
# arg1 = file, arg2 = file it depends on
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test ! -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# all package files
# only a few files have dependencies
for file in *.cpp *.h; do
if (test $file = "pair_cdeam.cpp") then
action pair_cdeam.cpp pair_eam_alloy.cpp
elif (test $file = "pair_cdeam.h") then
action pair_cdeam.h pair_eam_alloy.cpp
else
test -f ${file} && action $file
fi
done

1
src/USER-MISC/README
View File

@ -65,7 +65,6 @@ pair_style buck/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
pair_style coul/diel, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11
pair_style dipole/sf, Mario Orsi, orsimario at gmail.com, 8 Aug 11
pair_style edip, Luca Ferraro, luca.ferraro at caspur.it, 15 Sep 11
pair_style eam/cd, Alexander Stukowski, stukowski at mm.tu-darmstadt.de, 7 Nov 09
pair_style extep, Jaap Kroes (Radboud U), jaapkroes at gmail dot com, 28 Nov 17
pair_style gauss/cut, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11
pair_style lennard/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15

644
src/USER-MISC/pair_cdeam.cpp
View File

@ -1,644 +0,0 @@
/* ----------------------------------------------------------------------
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: Alexander Stukowski
Technical University of Darmstadt,
Germany Department of Materials Science
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "pair_cdeam.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
// This is for debugging purposes. The ASSERT() macro is used in the code to check
// if everything runs as expected. Change this to #if 0 if you don't need the checking.
#if 0
#define ASSERT(cond) ((!(cond)) ? my_failure(error,__FILE__,__LINE__) : my_noop())
inline void my_noop() {}
inline void my_failure(Error* error, const char* file, int line) {
char str[1024];
sprintf(str,"Assertion failure: File %s, line %i", file, line);
error->one(FLERR,str);
}
#else
#define ASSERT(cond)
#endif
#define MAXLINE 1024 // This sets the maximum line length in EAM input files.
PairCDEAM::PairCDEAM(LAMMPS *lmp, int _cdeamVersion) : PairEAM(lmp), PairEAMAlloy(lmp), cdeamVersion(_cdeamVersion)
{
single_enable = 0;
restartinfo = 0;
rhoB = NULL;
D_values = NULL;
hcoeff = NULL;
// Set communication buffer sizes needed by this pair style.
if(cdeamVersion == 1) {
comm_forward = 4;
comm_reverse = 3;
}
else if(cdeamVersion == 2) {
comm_forward = 3;
comm_reverse = 2;
}
else {
error->all(FLERR,"Invalid CD-EAM potential version.");
}
}
PairCDEAM::~PairCDEAM()
{
memory->destroy(rhoB);
memory->destroy(D_values);
if(hcoeff) delete[] hcoeff;
}
void PairCDEAM::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,rhoip,rhojp,recip,phi;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
// Grow per-atom arrays if necessary
if(atom->nmax > nmax) {
memory->destroy(rho);
memory->destroy(fp);
memory->destroy(rhoB);
memory->destroy(D_values);
nmax = atom->nmax;
memory->create(rho,nmax,"pair:rho");
memory->create(rhoB,nmax,"pair:rhoB");
memory->create(fp,nmax,"pair:fp");
memory->create(D_values,nmax,"pair:D_values");
}
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// Zero out per-atom arrays.
int m = nlocal + atom->nghost;
for(i = 0; i < m; i++) {
rho[i] = 0.0;
rhoB[i] = 0.0;
D_values[i] = 0.0;
}
// Stage I
// Compute rho and rhoB at each local atom site.
// Additionally calculate the D_i values here if we are using the one-site formulation.
// For the two-site formulation we have to calculate the D values in an extra loop (Stage II).
for(ii = 0; ii < inum; ii++) {
i = ilist[ii];
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];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if(rsq < cutforcesq) {
jtype = type[j];
double r = sqrt(rsq);
const EAMTableIndex index = radiusToTableIndex(r);
double localrho = RhoOfR(index, jtype, itype);
rho[i] += localrho;
if(jtype == speciesB) rhoB[i] += localrho;
if(newton_pair || j < nlocal) {
localrho = RhoOfR(index, itype, jtype);
rho[j] += localrho;
if(itype == speciesB) rhoB[j] += localrho;
}
if(cdeamVersion == 1 && itype != jtype) {
// Note: if the i-j interaction is not concentration dependent (because either
// i or j are not species A or B) then its contribution to D_i and D_j should
// be ignored.
// This if-clause is only required for a ternary.
if((itype == speciesA && jtype == speciesB) || (jtype == speciesA && itype == speciesB)) {
double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);
D_values[i] += Phi_AB;
if(newton_pair || j < nlocal)
D_values[j] += Phi_AB;
}
}
}
}
}
// Communicate and sum densities.
if(newton_pair) {
communicationStage = 1;
comm->reverse_comm_pair(this);
}
// fp = derivative of embedding energy at each atom
// phi = embedding energy at each atom
for(ii = 0; ii < inum; ii++) {
i = ilist[ii];
EAMTableIndex index = rhoToTableIndex(rho[i]);
fp[i] = FPrimeOfRho(index, type[i]);
if(eflag) {
phi = FofRho(index, type[i]);
if (eflag_global) eng_vdwl += phi;
if (eflag_atom) eatom[i] += phi;
}
}
// Communicate derivative of embedding function and densities
// and D_values (this for one-site formulation only).
communicationStage = 2;
comm->forward_comm_pair(this);
// The electron densities may not drop to zero because then the concentration would no longer be defined.
// But the concentration is not needed anyway if there is no interaction with another atom, which is the case
// if the electron density is exactly zero. That's why the following lines have been commented out.
//
//for(i = 0; i < nlocal + atom->nghost; i++) {
// if(rho[i] == 0 && (type[i] == speciesA || type[i] == speciesB))
// error->one(FLERR,"CD-EAM potential routine: Detected atom with zero electron density.");
//}
// Stage II
// This is only required for the original two-site formulation of the CD-EAM potential.
if(cdeamVersion == 2) {
// Compute intermediate value D_i for each atom.
for(ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
// This code line is required for ternary alloys.
if(itype != speciesA && itype != speciesB) continue;
double x_i = rhoB[i] / rho[i]; // Concentration at atom i.
for(jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
if(itype == jtype) continue;
// This code line is required for ternary alloys.
if(jtype != speciesA && jtype != speciesB) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if(rsq < cutforcesq) {
double r = sqrt(rsq);
const EAMTableIndex index = radiusToTableIndex(r);
// The concentration independent part of the cross pair potential.
double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);
// Average concentration of two sites
double x_ij = 0.5 * (x_i + rhoB[j]/rho[j]);
// Calculate derivative of h(x_ij) polynomial function.
double h_prime = evalHprime(x_ij);
D_values[i] += h_prime * Phi_AB / (2.0 * rho[i] * rho[i]);
if(newton_pair || j < nlocal)
D_values[j] += h_prime * Phi_AB / (2.0 * rho[j] * rho[j]);
}
}
}
// Communicate and sum D values.
if(newton_pair) {
communicationStage = 3;
comm->reverse_comm_pair(this);
}
communicationStage = 4;
comm->forward_comm_pair(this);
}
// Stage III
// Compute force acting on each atom.
for(ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
// Concentration at site i
double x_i = -1.0; // The value -1 indicates: no concentration dependence for all interactions of atom i.
// It will be replaced by the concentration at site i if atom i is either A or B.
double D_i, h_prime_i;
// This if-clause is only required for ternary alloys.
if((itype == speciesA || itype == speciesB) && rho[i] != 0.0) {
// Compute local concentration at site i.
x_i = rhoB[i]/rho[i];
ASSERT(x_i >= 0 && x_i<=1.0);
if(cdeamVersion == 1) {
// Calculate derivative of h(x_i) polynomial function.
h_prime_i = evalHprime(x_i);
D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]);
} else if(cdeamVersion == 2) {
D_i = D_values[i];
} else {
ASSERT(false);
}
}
for(jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if(rsq < cutforcesq) {
jtype = type[j];
double r = sqrt(rsq);
const EAMTableIndex index = radiusToTableIndex(r);
// rhoip = derivative of (density at atom j due to atom i)
// rhojp = derivative of (density at atom i due to atom j)
// psip needs both fp[i] and fp[j] terms since r_ij appears in two
// terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)
// hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
rhoip = RhoPrimeOfR(index, itype, jtype);
rhojp = RhoPrimeOfR(index, jtype, itype);
fpair = fp[i]*rhojp + fp[j]*rhoip;
recip = 1.0/r;
double x_j = -1; // The value -1 indicates: no concentration dependence for this i-j pair
// because atom j is not of species A nor B.
// This code line is required for ternary alloy.
if(jtype == speciesA || jtype == speciesB) {
ASSERT(rho[i] != 0.0);
ASSERT(rho[j] != 0.0);
// Compute local concentration at site j.
x_j = rhoB[j]/rho[j];
ASSERT(x_j >= 0 && x_j<=1.0);
double D_j=0.0;
if(cdeamVersion == 1) {
// Calculate derivative of h(x_j) polynomial function.
double h_prime_j = evalHprime(x_j);
D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]);
} else if(cdeamVersion == 2) {
D_j = D_values[j];
} else {
ASSERT(false);
}
double t2 = -rhoB[j];
if(itype == speciesB) t2 += rho[j];
fpair += D_j * rhoip * t2;
}
// This if-clause is only required for a ternary alloy.
// Actually we don't need it at all because D_i should be zero anyway if
// atom i has no concentration dependent interactions (because it is not species A or B).
if(x_i != -1.0) {
double t1 = -rhoB[i];
if(jtype == speciesB) t1 += rho[i];
fpair += D_i * rhojp * t1;
}
double phip;
double phi = PhiOfR(index, itype, jtype, recip, phip);
if(itype == jtype || x_i == -1.0 || x_j == -1.0) {
// Case of no concentration dependence.
fpair += phip;
} else {
// We have a concentration dependence for the i-j interaction.
double h=0.0;
if(cdeamVersion == 1) {
// Calculate h(x_i) polynomial function.
double h_i = evalH(x_i);
// Calculate h(x_j) polynomial function.
double h_j = evalH(x_j);
h = 0.5 * (h_i + h_j);
} else if(cdeamVersion == 2) {
// Average concentration.
double x_ij = 0.5 * (x_i + x_j);
// Calculate h(x_ij) polynomial function.
h = evalH(x_ij);
} else {
ASSERT(false);
}
fpair += h * phip;
phi *= h;
}
// Divide by r_ij and negate to get forces from gradient.
fpair /= -r;
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) evdwl = phi;
if(evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
}
}
}
if(vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairCDEAM::coeff(int narg, char **arg)
{
PairEAMAlloy::coeff(narg, arg);
// Make sure the EAM file is a CD-EAM binary alloy.
if(setfl->nelements < 2)
error->all(FLERR,"The EAM file must contain at least 2 elements to be used with the eam/cd pair style.");
// Read in the coefficients of the h polynomial from the end of the EAM file.
read_h_coeff(arg[2]);
// Determine which atom type is the A species and which is the B species in the alloy.
// By default take the first element (index 0) in the EAM file as the A species
// and the second element (index 1) in the EAM file as the B species.
speciesA = -1;
speciesB = -1;
for(int i = 1; i <= atom->ntypes; i++) {
if(map[i] == 0) {
if(speciesA >= 0)
error->all(FLERR,"The first element from the EAM file may only be mapped to a single atom type.");
speciesA = i;
}
if(map[i] == 1) {
if(speciesB >= 0)
error->all(FLERR,"The second element from the EAM file may only be mapped to a single atom type.");
speciesB = i;
}
}
if(speciesA < 0)
error->all(FLERR,"The first element from the EAM file must be mapped to exactly one atom type.");
if(speciesB < 0)
error->all(FLERR,"The second element from the EAM file must be mapped to exactly one atom type.");
}
/* ----------------------------------------------------------------------
Reads in the h(x) polynomial coefficients
------------------------------------------------------------------------- */
void PairCDEAM::read_h_coeff(char *filename)
{
if(comm->me == 0) {
// Open potential file
FILE *fptr;
char line[MAXLINE];
char nextline[MAXLINE];
fptr = force->open_potential(filename);
if (fptr == NULL) {
char str[128];
sprintf(str,"Cannot open EAM potential file %s", filename);
error->one(FLERR,str);
}
// h coefficients are stored at the end of the file.
// Skip to last line of file.
while(fgets(nextline, MAXLINE, fptr) != NULL) {
strcpy(line, nextline);
}
char* ptr = strtok(line, " \t\n\r\f");
int degree = atoi(ptr);
nhcoeff = degree+1;
hcoeff = new double[nhcoeff];
int i = 0;
while((ptr = strtok(NULL," \t\n\r\f")) != NULL && i < nhcoeff) {
hcoeff[i++] = atof(ptr);
}
if(i != nhcoeff || nhcoeff < 1)
error->one(FLERR,"Failed to read h(x) function coefficients from EAM file.");
// Close the potential file.
fclose(fptr);
}
MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world);
if(comm->me != 0) hcoeff = new double[nhcoeff];
MPI_Bcast(hcoeff, nhcoeff, MPI_DOUBLE, 0, world);
}
/* ---------------------------------------------------------------------- */
int PairCDEAM::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
m = 0;
if(communicationStage == 2) {
if(cdeamVersion == 1) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = fp[j];
buf[m++] = rho[j];
buf[m++] = rhoB[j];
buf[m++] = D_values[j];
}
return m;
}
else if(cdeamVersion == 2) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = fp[j];
buf[m++] = rho[j];
buf[m++] = rhoB[j];
}
return m;
}
else { ASSERT(false); return 0; }
}
else if(communicationStage == 4) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = D_values[j];
}
return m;
}
else return 0;
}
/* ---------------------------------------------------------------------- */
void PairCDEAM::unpack_forward_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if(communicationStage == 2) {
if(cdeamVersion == 1) {
for(i = first; i < last; i++) {
fp[i] = buf[m++];
rho[i] = buf[m++];
rhoB[i] = buf[m++];
D_values[i] = buf[m++];
}
}
else if(cdeamVersion == 2) {
for(i = first; i < last; i++) {
fp[i] = buf[m++];
rho[i] = buf[m++];
rhoB[i] = buf[m++];
}
} else {
ASSERT(false);
}
}
else if(communicationStage == 4) {
for(i = first; i < last; i++) {
D_values[i] = buf[m++];
}
}
}
/* ---------------------------------------------------------------------- */
int PairCDEAM::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if(communicationStage == 1) {
if(cdeamVersion == 1) {
for(i = first; i < last; i++) {
buf[m++] = rho[i];
buf[m++] = rhoB[i];
buf[m++] = D_values[i];
}
return m;
}
else if(cdeamVersion == 2) {
for(i = first; i < last; i++) {
buf[m++] = rho[i];
buf[m++] = rhoB[i];
}
return m;
}
else { ASSERT(false); return 0; }
}
else if(communicationStage == 3) {
for(i = first; i < last; i++) {
buf[m++] = D_values[i];
}
return m;
}
else return 0;
}
/* ---------------------------------------------------------------------- */
void PairCDEAM::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
if(communicationStage == 1) {
if(cdeamVersion == 1) {
for(i = 0; i < n; i++) {
j = list[i];
rho[j] += buf[m++];
rhoB[j] += buf[m++];
D_values[j] += buf[m++];
}
} else if(cdeamVersion == 2) {
for(i = 0; i < n; i++) {
j = list[i];
rho[j] += buf[m++];
rhoB[j] += buf[m++];
}
} else {
ASSERT(false);
}
}
else if(communicationStage == 3) {
for(i = 0; i < n; i++) {
j = list[i];
D_values[j] += buf[m++];
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairCDEAM::memory_usage()
{
double bytes = 2 * nmax * sizeof(double);
return PairEAMAlloy::memory_usage() + bytes;
}

36
src/fix_nh.cpp
View File

@ -798,7 +798,7 @@ void FixNH::setup(int vflag)
if (pstat_flag) {
double kt = boltz * t_target;
double nkt = atom->natoms * kt;
double nkt = (atom->natoms + 1) * kt;
for (int i = 0; i < 3; i++)
if (p_flag[i])
@ -1446,7 +1446,7 @@ double FixNH::compute_scalar()
double volume;
double energy;
double kt = boltz * t_target;
double lkt_press = kt;
double lkt_press = 0.0;
int ich;
if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
else volume = domain->xprd * domain->yprd;
@ -1477,15 +1477,21 @@ double FixNH::compute_scalar()
// sum is over barostatted dimensions
if (pstat_flag) {
for (i = 0; i < 3; i++)
if (p_flag[i])
for (i = 0; i < 3; i++) {
if (p_flag[i]) {
energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i] +
p_hydro*(volume-vol0) / (pdim*nktv2p);
lkt_press += kt;
}
}
if (pstyle == TRICLINIC) {
for (i = 3; i < 6; i++)
if (p_flag[i])
for (i = 3; i < 6; i++) {
if (p_flag[i]) {
energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i];
lkt_press += kt;
}
}
}
// extra contributions from thermostat chain for barostat
@ -1818,15 +1824,15 @@ void FixNH::nhc_temp_integrate()
void FixNH::nhc_press_integrate()
{
int ich,i;
int ich,i,pdof;
double expfac,factor_etap,kecurrent;
double kt = boltz * t_target;
double lkt_press = kt;
double lkt_press;
// Update masses, to preserve initial freq, if flag set
if (omega_mass_flag) {
double nkt = atom->natoms * kt;
double nkt = (atom->natoms + 1) * kt;
for (int i = 0; i < 3; i++)
if (p_flag[i])
omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
@ -1850,14 +1856,22 @@ void FixNH::nhc_press_integrate()
}
kecurrent = 0.0;
pdof = 0;
for (i = 0; i < 3; i++)
if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
if (p_flag[i]) {
kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
pdof++;
}
if (pstyle == TRICLINIC) {
for (i = 3; i < 6; i++)
if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
if (p_flag[i]) {
kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
pdof++;
}
}
lkt_press = pdof * kt;
etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
double ncfac = 1.0/nc_pchain;

2
src/version.h
View File

@ -1 +1 @@
#define LAMMPS_VERSION "16 Jul 2018"
#define LAMMPS_VERSION "2 Aug 2018"

82
tools/emacs/README.md Normal file
View File

@ -0,0 +1,82 @@
# GNU Emacs Syntax Highlighting
> Copyright (C) 2010-2018 Aidan Thompson <athomps at sandia.gov>
> Copyright (C) 2018 Rohit Goswami <r95g10 at gmail.com>
The `lammps-mode.el` file provided in this directory will enable syntax
highlighting for the lammps script syntax in GNU Emacs. The groupings of
commands were originally copied from `tools/vim`.
## Installation
**Requirements: GNU Emacs 24.\***
### Obtaining the Package
#### MELPA
The easiest installation method is via MELPA and it is advisable to use one of
the many [MELPA installation methods](https://melpa.org/#/getting-started).
For example, with [use-package](https://github.com/jwiegley/use-package) one can
simply use the following:
``` emacs-lisp
(use-package lammps-mode)
```
#### Manually
Assuming for some reason you have downloaded the file to `~/.emacs.d/lisp` you
would do the following (kanged [from here](http://ergoemacs.org/emacs/emacs_installing_packages.html)):
``` emacs-lisp
;; Tell emacs where is your personal elisp lib dir
(add-to-list 'load-path "~/.emacs.d/lisp/")
;; load the package.
(load "lammps-mode")
```
### Autoloading \& Auto-recognition
To automatically turn on the LAMMPS mode for editing your input scripts,
use the following line as the **first** line of your script:
```
# -*- lammps -*-
```
For automatically switching on the LAMMPS mode based on filename patterns,
e.g. for `in.*` and `*.lmp` files, add the following code to your `.emacs`:
``` emacs-lisp
(autoload 'lammps-mode "lammps-mode.el" "LAMMPS mode." t)
(setq auto-mode-alist (append auto-mode-alist
'(("in\\." . lammps-mode))
'(("\\.lmp\\'" . lammps-mode))
))
```
## Status
By far not all commands are included in the syntax file (lammps-mode.el). You
can easily add new ones to the existing classes.
## Implementation Details
`lammps-mode` is derived from `shell-script-mode` which provides some basic
syntax highlighting of strings, comments, etc.
The MELPA recipe used for this package is simply:
``` emacs-lisp
(lammps-mode :fetcher github :repo "HaoZeke/lammps-mode")
```
## Caveats
* Does not work with Xemacs [See [this comment](https://github.com/lammps/lammps/pull/1022#issuecomment-408871233)]
## License
[GNU GPL v2](https://github.com/HaoZeke/lammps-mode/blob/master/LICENSE).
Check the file for more details.

23
tools/emacs/README.txt
View File

@ -1,23 +0,0 @@
=== Emacs Syntax Highlighting ===
Created by Aidan Thompson 12/2010
===============================
The lammps.el file provided in this directory will enable syntax
highlighting for the lammps script syntax in emacs. The groupings
of commands were copied from tools/vim. The simulation scripts have to
end on *.lmp or start with in.* (see lammps.el). By far not all
commands are included in the syntax file (lammps.el).
You can easily add new ones to the existing classes.
'lammps-mode' is derived from 'shell-script-mode' which provides
some basic syntax highlighting of strings, comments, etc.
=To enable the highlighting:
============================
(0) Create/edit the emacs init file ~/.emacs to contain:
(load "~/.emacs.d/lammps")
This file may also be called ~/.emacs.el, or ~/.emacs.d/init.el
(1) Copy lammps.el to the directory ~/.emacs.d

59
tools/emacs/lammps.el → tools/emacs/lammps-mode.el
View File

@ -1,7 +1,48 @@
;; LAMMPS auto-mode
;;; lammps-mode.el --- basic syntax highlighting for LAMMPS files
;; Copyright (C) 2010-18 Aidan Thompson
;; Copyright (C) 2018 Rohit Goswami
;; Author: Aidan Thompson <athomps at sandia.gov>
;; Maintainer: Rohit Goswami <r95g10 at gmail.com>
;; Created: December 4, 2010
;; Modified: July 30, 2018
;; Version: 1.5.0
;; Keywords: languages, faces
;; Homepage: https://github.com/lammps/lammps/tree/master/tools/emacs
;; Package-Requires: ((emacs "24.4"))
;; This file is not part of GNU Emacs.
;; This program is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 2 of the License, or
;; (at your option) any later version.
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
;; GNU General Public License for more details.
;; You should have received a copy of the GNU General Public License along
;; with this program; if not, write to the Free Software Foundation, Inc.,
;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
;;; Commentary:
;; translation of keyword classes from tools/vim
;; see http://xahlee.org/emacs/elisp_syntax_coloring.html
;; Put this in your .emacs file to enable autoloading of lammps-mode
;; and auto-recognition of "in.*" and "*.lmp" files:
;;
;; (autoload 'lammps-mode "lammps-mode.el" "LAMMPS mode." t)
;; (setq auto-mode-alist (append auto-mode-alist
;; '(("in\\." . lammps-mode))
;; '(("\\.lmp\\'" . lammps-mode))
;; ))
;;
;;; Code:
;; define several keyword classes
(defvar lammps-output
'("log"
@ -136,6 +177,8 @@
(defvar lammps-variable-regexp
"\\$\\({[a-zA-Z0-9_]+}\\)\\|\\$[A-Za-z]")
(defvar lammps-font-lock-keywords)
;; clear memory
(setq lammps-output nil)
(setq lammps-read nil)
@ -151,8 +194,7 @@
;; create the list for font-lock.
;; each class of keyword is given a particular face
(setq
lammps-font-lock-keywords
(setq lammps-font-lock-keywords
`((,lammps-output-regexp . font-lock-function-name-face)
(,lammps-read-regexp . font-lock-preprocessor-face)
(,lammps-lattice-regexp . font-lock-type-face)
@ -199,12 +241,5 @@
(setq lammps-comment-regexp nil)
(setq lammps-variable-regexp nil))
;; apply it to specified filename patterns
(setq
auto-mode-alist
(append
auto-mode-alist
'(("in\\." . lammps-mode))
'(("\\.lmp\\'" . lammps-mode))
))
(provide 'lammps-mode)
;;; lammps-mode.el ends here