lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@939 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2007-10-03 16:28:14 +00:00
parent 10cb6b72e9
commit c65b2445eb
19 changed files with 4360 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
src/USER-ACKLAND/Install.csh Normal file
View File

@ -0,0 +1,20 @@
# Install/unInstall package classes in LAMMPS
if ($1 == 1) then
cp -p style_user_ackland.h ..
cp -p compute_ackland_atom.cpp ..
cp -p compute_ackland_atom.h ..
else if ($1 == 0) then
rm ../style_user_ackland.h
touch ../style_user_ackland.h
rm ../compute_ackland_atom.cpp
rm ../compute_ackland_atom.h
endif

18
src/USER-ACKLAND/README Normal file
View File

@ -0,0 +1,18 @@
The files in this directory are a user-contributed package for LAMMPS.
The person who created these files is Gerolf Ziegenhain
(gerolf@ziegenhain.com). Contact him directly if you have questions.
This package implements a "compute ackland/atom" command which can be
used in a LAMMPS input script. Like other per-atom compute commands,
the results can be accessed when dumping atom information to a file,
or by other fixes that do averaging of various kinds. See the
documentation files for these commands for details.
The Ackland computation is a means of detecting local lattice
structure around an atom, as described in G. Ackland,
PRB(2006)73:054104.
The output is a number with the following mapping:
enum{UNKNOWN,BCC,FCC,HCP,ICO};

385
src/USER-ACKLAND/compute_ackland_atom.cpp Normal file
View File

@ -0,0 +1,385 @@
/* ----------------------------------------------------------------------
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: G. Ziegenhain, gerolf@ziegenhain.com
Copyright (C) 2007
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_ackland_atom.h"
#include "atom.h"
#include "modify.h"
#include "update.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "pair.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include <math.h>
using namespace LAMMPS_NS;
enum{UNKNOWN,BCC,FCC,HCP,ICO};
/* ---------------------------------------------------------------------- */
ComputeAcklandAtom::ComputeAcklandAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all("Illegal compute ackland/atom command");
peratom_flag = 1;
size_peratom = 0;
nmax = 0;
structure = NULL;
maxneigh = 0;
distsq = NULL;
nearest = NULL;
nearest_n0 = NULL;
nearest_n1 = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeAcklandAtom::~ComputeAcklandAtom()
{
memory->sfree(structure);
memory->sfree(distsq);
memory->sfree(nearest);
memory->sfree(nearest_n0);
memory->sfree(nearest_n1);
}
/* ---------------------------------------------------------------------- */
void ComputeAcklandAtom::init()
{
// need an occasional full neighbor list
int irequest = neighbor->request((void *) this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->occasional = 1;
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"ackland/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning("More than one compute ackland/atom");
}
/* ---------------------------------------------------------------------- */
void ComputeAcklandAtom::compute_peratom()
{
int i,j,ii,jj,k,n,inum,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,value;
int *ilist,*jlist,*numneigh,**firstneigh;
double pairs[66];
int chi[8];
// grow structure array if necessary
if (atom->nlocal > nmax) {
memory->sfree(structure);
nmax = atom->nmax;
structure = (double *)
memory->smalloc(nmax*sizeof(double),"compute/ackland/atom:ackland");
scalar_atom = structure;
}
// invoke half neighbor list (will copy or build if necessary)
neighbor->build_one(list->index);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// compute structure parameter for each atom in group
// use full neighbor list
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;
double cutsq = force->pair->cutforce * force->pair->cutforce;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
// ensure distsq and nearest arrays are long enough
if (jnum > maxneigh) {
memory->sfree(distsq);
memory->sfree(nearest);
memory->sfree(nearest_n0);
memory->sfree(nearest_n1);
maxneigh = jnum;
distsq = (double *) memory->smalloc(maxneigh*sizeof(double),
"compute/ackland/atom:distsq");
nearest = (int *) memory->smalloc(maxneigh*sizeof(int),
"compute/ackland/atom:nearest");
nearest_n0 = (int *) memory->smalloc(maxneigh*sizeof(int),
"compute/ackland/atom:nearest_n0");
nearest_n1 = (int *) memory->smalloc(maxneigh*sizeof(int),
"compute/ackland/atom:nearest_n1");
}
// loop over list of all neighbors within force cutoff
// distsq[] = distance sq to each
// nearest[] = atom indices of neighbors
n = 0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
if (j >= nall) j %= nall;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutsq) {
distsq[n] = rsq;
nearest[n++] = j;
}
}
// Select 6 nearest neighbors
select2(6,n,distsq,nearest);
// Mean squared separation
double r0_sq = 0.;
for (j = 0; j < 6; j++)
r0_sq += distsq[j];
r0_sq /= 6.;
// n0 near neighbors with: distsq<1.45*r0_sq
// n1 near neighbors with: distsq<1.55*r0_sq
double n0_dist_sq = 1.45*r0_sq,
n1_dist_sq = 1.55*r0_sq;
int n0 = 0, n1 = 0;
for (j = 0; j < n; j++) {
if (distsq[j] < n1_dist_sq) {
nearest_n1[n1++] = nearest[j];
if (distsq[j] < n0_dist_sq) {
nearest_n0[n0++] = nearest[j];
}
}
}
// Evaluate all angles <(r_ij,rik) forall n0 particles with: distsq<1.45*r0_sq
double bond_angle;
double norm_j, norm_k;
chi[0] = chi[1] = chi[2] = chi[3] = chi[4] = chi[5] = chi[6] = chi[7] = 0;
double x_ij, y_ij, z_ij, x_ik, y_ik, z_ik;
for (j = 0; j < n0; j++) {
x_ij = x[i][0]-x[nearest_n0[j]][0];
y_ij = x[i][1]-x[nearest_n0[j]][1];
z_ij = x[i][2]-x[nearest_n0[j]][2];
norm_j = sqrt (x_ij*x_ij + y_ij*y_ij + z_ij*z_ij);
if (norm_j <= 0.) continue;
for (k = j+1; k < n0; k++) {
x_ik = x[i][0]-x[nearest_n0[k]][0];
y_ik = x[i][1]-x[nearest_n0[k]][1];
z_ik = x[i][2]-x[nearest_n0[k]][2];
norm_k = sqrt (x_ik*x_ik + y_ik*y_ik + z_ik*z_ik);
if (norm_k <= 0.)
continue;
bond_angle = (x_ij*x_ik + y_ij*y_ik + z_ij*z_ik) / (norm_j*norm_k);
// Histogram for identifying the relevant peaks
if (-1. <= bond_angle && bond_angle < -0.945) { chi[0]++; }
else if (-0.945 <= bond_angle && bond_angle < -0.915) { chi[1]++; }
else if (-0.915 <= bond_angle && bond_angle < -0.755) { chi[2]++; }
else if (-0.755 <= bond_angle && bond_angle < -0.195) { chi[3]++; }
else if (-0.195 <= bond_angle && bond_angle < 0.195) { chi[4]++; }
else if (0.195 <= bond_angle && bond_angle < 0.245) { chi[5]++; }
else if (0.245 <= bond_angle && bond_angle < 0.795) { chi[6]++; }
else if (0.795 <= bond_angle && bond_angle < 1.) { chi[7]++; }
}
}
// Deviations from the different lattice structures
double delta_bcc = 0.35*chi[4]/(double)(chi[5]+chi[6]-chi[4]),
delta_cp = fabs(1.-(double)chi[6]/24.),
delta_fcc = 0.61*(fabs((double)(chi[0]+chi[1]-6.))+(double)chi[2])/6.,
delta_hcp = (fabs((double)chi[0]-3.)+fabs((double)chi[0]+(double)chi[1]+(double)chi[2]+(double)chi[3]-9.))/12.;
// Identification of the local structure according to the reference
if (chi[0] == 7) { delta_bcc = 0.; }
else if (chi[0] == 6) { delta_fcc = 0.; }
else if (chi[0] <= 3) { delta_hcp = 0.; }
if (chi[7] > 0.)
structure[i] = UNKNOWN;
else
if (chi[4] < 3.)
{
if (n1 > 13 || n1 < 11)
structure[i] = UNKNOWN;
else
structure[i] = ICO;
} else
if (delta_bcc <= delta_cp)
{
if (n1 < 11)
structure[i] = UNKNOWN;
else
structure[i] = BCC;
} else
if (n1 > 12 || n1 < 11)
structure[i] = UNKNOWN;
else
if (delta_fcc < delta_hcp)
structure[i] = FCC;
else
structure[i] = HCP;
} // end loop over all particles
}
}
/* ----------------------------------------------------------------------
2 select routines from Numerical Recipes (slightly modified)
find k smallest values in array of length n
2nd routine sorts auxiliary array at same time
------------------------------------------------------------------------- */
#define SWAP(a,b) tmp = a; a = b; b = tmp;
#define ISWAP(a,b) itmp = a; a = b; b = itmp;
void ComputeAcklandAtom::select(int k, int n, double *arr)
{
int i,ir,j,l,mid;
double a,tmp;
arr--;
l = 1;
ir = n;
for (;;) {
if (ir <= l+1) {
if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir])
}
return;
} else {
mid=(l+ir) >> 1;
SWAP(arr[mid],arr[l+1])
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir])
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir])
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1])
}
i = l+1;
j = ir;
a = arr[l+1];
for (;;) {
do i++; while (arr[i] < a);
do j--; while (arr[j] > a);
if (j < i) break;
SWAP(arr[i],arr[j])
}
arr[l+1] = arr[j];
arr[j] = a;
if (j >= k) ir = j-1;
if (j <= k) l = i;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputeAcklandAtom::select2(int k, int n, double *arr, int *iarr)
{
int i,ir,j,l,mid,ia,itmp;
double a,tmp;
arr--;
iarr--;
l = 1;
ir = n;
for (;;) {
if (ir <= l+1) {
if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir])
ISWAP(iarr[l],iarr[ir])
}
return;
} else {
mid=(l+ir) >> 1;
SWAP(arr[mid],arr[l+1])
ISWAP(iarr[mid],iarr[l+1])
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir])
ISWAP(iarr[l],iarr[ir])
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir])
ISWAP(iarr[l+1],iarr[ir])
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1])
ISWAP(iarr[l],iarr[l+1])
}
i = l+1;
j = ir;
a = arr[l+1];
ia = iarr[l+1];
for (;;) {
do i++; while (arr[i] < a);
do j--; while (arr[j] > a);
if (j < i) break;
SWAP(arr[i],arr[j])
ISWAP(iarr[i],iarr[j])
}
arr[l+1] = arr[j];
arr[j] = a;
iarr[l+1] = iarr[j];
iarr[j] = ia;
if (j >= k) ir = j-1;
if (j <= k) l = i;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
int ComputeAcklandAtom::memory_usage()
{
int bytes = nmax * sizeof(double);
return bytes;
}

42
src/USER-ACKLAND/compute_ackland_atom.h Normal file
View File

@ -0,0 +1,42 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifndef COMPUTE_ACKLAND_ATOM_H
#define COMPUTE_ACKLAND_ATOM_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeAcklandAtom : public Compute {
public:
ComputeAcklandAtom(class LAMMPS *, int, char **);
~ComputeAcklandAtom();
void init();
void compute_peratom();
int memory_usage();
private:
int nmax,maxneigh;
double *distsq;
int *nearest, *nearest_n0, *nearest_n1;
double *structure;
class NeighList *list;
void select(int, int, double *);
void select2(int, int, double *, int *);
};
}
#endif

20
src/USER-ACKLAND/style_user_ackland.h Normal file
View File

@ -0,0 +1,20 @@
/* ----------------------------------------------------------------------
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 ComputeInclude
#include "compute_ackland_atom.h"
#endif
#ifdef ComputeClass
ComputeStyle(ackland/atom,ComputeAcklandAtom)
#endif

34
src/USER-EWALDN/Install.csh Normal file
View File

@ -0,0 +1,34 @@
# Install/unInstall package classes in LAMMPS
if ($1 == 1) then
cp -p style_user_ewaldn.h ..
cp -p ewald_n.cpp ..
cp -p pair_buck_coul.cpp ..
cp -p pair_lj_coul.cpp ..
cp -p ewald_n.h ..
cp -p pair_buck_coul.h ..
cp -p pair_lj_coul.h ..
cp -p math_vector.h ..
cp -p math_complex.h ..
else if ($1 == 0) then
rm ../style_user_ewaldn.h
touch ../style_user_ewaldn.h
rm ../ewald_n.cpp
rm ../pair_buck_coul.cpp
rm ../pair_lj_coul.cpp
rm ../ewald_n.h
rm ../pair_buck_coul.h
rm ../pair_lj_coul.h
rm ../math_vector.h
rm ../math_complex.h
endif

22
src/USER-EWALDN/README Normal file
View File

@ -0,0 +1,22 @@
The files in this directory are a user-contributed package for LAMMPS.
The person who created these files is Pieter in' t Veld at Sandia
(pjintve@sandia.gov). Contact him directly if you have questions.
This package implements 3 commands which can be used in a LAMMPS input
script: pair_style lj/coul, pair_style buck/coul, and kspace_style
ewald/n. See the documentation files for these commands for details.
The "kspace_style ewald/n" command is similar to standard Ewald for
charges, but also enables the Lennard-Jones interaction, or any 1/r^N
interaction to be of infinite extent, instead of being cutoff. LAMMPS
pair potentials for long-range Coulombic interactions, such as
lj/cut/coul/long can be used with ewald/n. The two new pair_style
commands provide the modifications for the short-range LJ and
Buckingham interactions that can also be used with ewald/n.
Another advantage of kspace_style ewald/n is that it can be used with
non-orthogonal (triclinic symmetry) simulation boxes, either for just
long-range Coulombic interactions, or for both Coulombic and 1/r^N LJ
or Buckingham, which is not currently possible for other kspace styles
such as PPPM and ewald.

704
src/USER-EWALDN/ewald_n.cpp Normal file
View File

@ -0,0 +1,704 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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: Pieter J. in 't Veld (SNL)
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
module: ewald_n.cpp
author: Pieter J. in 't Veld for SNL
date: September 13, 2006.
usage: kspace ewald/n precision
precision set precision of all orders
remarks: - cut off and precision are assumed identical for all orders
- coulombics from Macromolecules 1989, 93, 7320
- dipoles from J. Chem. Phys. 2000, 113, 10913
* ---------------------------------------------------------------------- */
#include "mpi.h"
#include "string.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "ewald_n.h"
#include "math_vector.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define KSPACE_ILLEGAL "Illegal kspace_style ewald/n command"
#define KSPACE_ORDER "Unsupported order in kspace_style ewald/n for"
#define KSPACE_MIX "Unsupported mixing rule in kspace_style ewald/n for"
enum{GEOMETRIC,ARITHMETIC,SIXTHPOWER}; // same as in pair.cpp
//#define DEBUG
/* ---------------------------------------------------------------------- */
EwaldN::EwaldN(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
{
if (narg!=1) error->all(KSPACE_ILLEGAL);
precision = fabs(atof(arg[0]));
memset(function, 0, EWALD_NORDER*sizeof(int));
kenergy = kvirial = NULL;
cek_local = cek_global = NULL;
ekr_local = NULL;
hvec = NULL;
kvec = NULL;
B = NULL;
first_output = 0;
}
EwaldN::~EwaldN()
{
deallocate();
delete [] ekr_local;
delete [] B;
}
/* --------------------------------------------------------------------- */
void EwaldN::init()
{
nkvec = nkvec_max = nevec = nevec_max = bytes = 0;
nfunctions = nsums = sums = 0;
nbox = -1;
if (!comm->me) { // output message
if (screen) fprintf(screen,"EwaldN initialization ...\n");
if (logfile) fprintf(logfile,"EwaldN initialization ...\n");
}
if (domain->dimension == 2) // check for errors
error->all("Cannot use EwaldN with 2d simulation");
if (slabflag == 0 && domain->nonperiodic > 0)
error->all("Cannot use nonperiodic boundaries with EwaldN");
if (slabflag == 1) {
if (domain->xperiodic != 1 || domain->yperiodic != 1 ||
domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1)
error->all("Incorrect boundaries with slab EwaldN");
}
qqrd2e = force->qqrd2e; // check pair_style
//mumurd2e = force->mumurd2e;
//dielectric = force->dielectric;
mumurd2e = dielectric = 1.0;
Pair *pair = force->pair;
void *ptr = pair ? pair->extract_ptr("ewald_order") : NULL;
if (!ptr) error->all("KSpace style is incompatible with Pair style");
int ewald_order = *((int *) ptr);
int ewald_mix = *((int *) pair->extract_ptr("ewald_mix"));
double cutoff = *((double *) pair->extract_ptr("ewald_cut"));
memset(function, 0, EWALD_NFUNCS*sizeof(int));
for (int i=0; i<=EWALD_NORDER; ++i) // transcribe order
if (ewald_order&(1<<i)) { // from pair_style
int n[] = EWALD_NSUMS, k;
char str[128];
switch (i) {
case 1:
k = 0; break;
case 3:
k = 3; break;
case 6:
if (ewald_mix==GEOMETRIC) { k = 1; break; }
else if (ewald_mix==ARITHMETIC) { k = 2; break; }
sprintf(str, "%s pair_style %s", KSPACE_MIX, force->pair_style);
error->all(str);
default:
sprintf(str, "%s pair_style %s", KSPACE_ORDER, force->pair_style);
error->all(str);
}
nfunctions += function[k] = 1;
nsums += n[k];
}
g_ewald = (1.35 - 0.15*log(precision))/cutoff; // determine resolution
g2_max = -4.0*g_ewald*g_ewald*log(precision);
if (!comm->me) { // output results
if (screen) fprintf(screen, " G vector = %g\n", g_ewald);
if (logfile) fprintf(logfile, " G vector = %g\n", g_ewald);
}
}
/* ----------------------------------------------------------------------
adjust EwaldN coeffs, called initially and whenever volume has changed
------------------------------------------------------------------------- */
void EwaldN::setup()
{
volume = shape_det(domain->h)*slab_volfactor; // cell volume
memcpy(unit, domain->h_inv, sizeof(shape)); // wave vector units
shape_scalar_mult(unit, 2.0*M_PI);
unit[2] /= slab_volfactor;
//int nbox_old = nbox, nkvec_old = nkvec;
if (precision>=1) nbox = 0;
else {
vector n = {1.0, 1.0, 1.0}; // based on cutoff
vec_scalar_mult(n, g_ewald*sqrt(-log(precision))/M_PI);
shape_vec_dot(n, n, domain->h);
n[2] *= slab_volfactor;
nbox = (int) n[0];
if (nbox<(int) n[1]) nbox = (int) n[1];
if (nbox<(int) n[2]) nbox = (int) n[2];
}
reallocate();
coefficients(); // compute coeffs
init_coeffs();
init_coeff_sums();
init_self();
if (!(first_output||comm->me)) { // output on first
first_output = 1;
if (screen) fprintf(screen,
" vectors: nbox = %d, nkvec = %d\n", nbox, nkvec);
if (logfile) fprintf(logfile,
" vectors: nbox = %d, nkvec = %d\n", nbox, nkvec);
}
}
void EwaldN::reallocate() // allocate memory
{
int ix, iy, iz;
int nkvec_max = nkvec;
vector h;
nkvec = 0; // determine size(kvec)
int kflag[(nbox+1)*(2*nbox+1)*(2*nbox+1)], *flag = kflag;
for (ix=0; ix<=nbox; ++ix)
for (iy=-nbox; iy<=nbox; ++iy)
for (iz=-nbox; iz<=nbox; ++iz)
if (!(ix||iy||iz)) *(flag++) = 0;
else if ((!ix)&&(iy<0)) *(flag++) = 0;
else if ((!(ix||iy))&&(iz<0)) *(flag++) = 0; // use symmetry
else {
h[0] = unit[0]*ix;
h[1] = unit[5]*ix+unit[1]*iy;
h[2] = unit[4]*ix+unit[3]*iy+unit[2]*iz;
if ((*(flag++) = h[0]*h[0]+h[1]*h[1]+h[2]*h[2]<=g2_max)) ++nkvec;
}
if (nkvec>nkvec_max) {
deallocate(); // free memory
hvec = new hvector[nkvec]; // hvec
bytes += (nkvec-nkvec_max)*sizeof(hvector);
kvec = new kvector[nkvec]; // kvec
bytes += (nkvec-nkvec_max)*sizeof(kvector);
kenergy = new double[nkvec*nfunctions]; // kenergy
bytes += (nkvec-nkvec_max)*nfunctions*sizeof(double);
kvirial = new double[6*nkvec*nfunctions]; // kvirial
bytes += 6*(nkvec-nkvec_max)*nfunctions*sizeof(double);
cek_local = new complex[nkvec*nsums]; // cek_local
bytes += (nkvec-nkvec_max)*nsums*sizeof(complex);
cek_global = new complex[nkvec*nsums]; // cek_global
bytes += (nkvec-nkvec_max)*nsums*sizeof(complex);
nkvec_max = nkvec;
}
flag = kflag; // create index and
kvector *k = kvec; // wave vectors
hvector *hi = hvec;
for (ix=0; ix<=nbox; ++ix)
for (iy=-nbox; iy<=nbox; ++iy)
for (iz=-nbox; iz<=nbox; ++iz)
if (*(flag++)) {
hi->x = unit[0]*ix;
hi->y = unit[5]*ix+unit[1]*iy;
(hi++)->z = unit[4]*ix+unit[3]*iy+unit[2]*iz;
k->x = ix+nbox; k->y = iy+nbox; (k++)->z = iz+nbox; }
}
void EwaldN::reallocate_atoms()
{
if ((nevec = atom->nmax*(2*nbox+1))<=nevec_max) return;
delete [] ekr_local;
ekr_local = new cvector[nevec];
bytes += (nevec-nevec_max)*sizeof(cvector);
nevec_max = nevec;
}
void EwaldN::deallocate() // free memory
{
delete [] hvec; hvec = NULL;
delete [] kvec; kvec = NULL;
delete [] kenergy; kenergy = NULL;
delete [] kvirial; kvirial = NULL;
delete [] cek_local; cek_local = NULL;
delete [] cek_global; cek_global = NULL;
}
void EwaldN::coefficients() // set up pre-factors
{
vector h;
hvector *hi = hvec, *nh;
double eta2 = 0.25/(g_ewald*g_ewald);
double b1, b2, expb2, h1, h2, c1, c2;
double *ke = kenergy, *kv = kvirial;
int func0 = function[0], func12 = function[1]||function[2],
func3 = function[3];
for (nh = (hi = hvec)+nkvec; hi<nh; ++hi) { // wave vectors
memcpy(h, hi, sizeof(vector));
expb2 = exp(-(b2 = (h2 = vec_dot(h, h))*eta2));
if (func0) { // qi*qj/r coeffs
*(ke++) = c1 = expb2/h2;
*(kv++) = c1-(c2 = 2.0*c1*(1.0+b2)/h2)*h[0]*h[0];
*(kv++) = c1-c2*h[1]*h[1]; // lammps convention
*(kv++) = c1-c2*h[2]*h[2]; // instead of voigt
*(kv++) = -c2*h[1]*h[0];
*(kv++) = -c2*h[2]*h[0];
*(kv++) = -c2*h[2]*h[1];
}
if (func12) { // -Bij/r^6 coeffs
b1 = sqrt(b2); // minus sign folded
h1 = sqrt(h2); // into constants
*(ke++) = c1 = -h1*h2*((c2=sqrt(M_PI)*erfc(b1))+(0.5/b2-1.0)*expb2/b1);
*(kv++) = c1-(c2 = 3.0*h1*(c2-expb2/b1))*h[0]*h[0];
*(kv++) = c1-c2*h[1]*h[1]; // lammps convention
*(kv++) = c1-c2*h[2]*h[2]; // instead of voigt
*(kv++) = -c2*h[1]*h[0];
*(kv++) = -c2*h[2]*h[0];
*(kv++) = -c2*h[2]*h[1];
}
if (func3) { // dipole coeffs
*(ke++) = c1 = expb2/h2;
*(kv++) = c1-(c2 = 2.0*c1*(1.0+b2)/h2)*h[0]*h[0];
*(kv++) = c1-c2*h[1]*h[1]; // lammps convention
*(kv++) = c1-c2*h[2]*h[2]; // instead of voigt
*(kv++) = -c2*h[1]*h[0];
*(kv++) = -c2*h[2]*h[0];
*(kv++) = -c2*h[2]*h[1];
}
}
}
void EwaldN::init_coeffs() // local pair coeffs
{
int n = atom->ntypes;
if (function[1]) { // geometric 1/r^6
double **b = (double **) force->pair->extract_ptr("B");
delete [] B;
B = new double[n+1];
bytes += (n+1)*sizeof(double);
for (int i=0; i<=n; ++i) B[i] = sqrt(fabs(b[i][i]));
}
if (function[2]) { // arithmetic 1/r^6
double **epsilon = (double **) force->pair->extract_ptr("epsilon");
double **sigma = (double **) force->pair->extract_ptr("sigma");
if (!(epsilon&&sigma))
error->all("epsilon or sigma reference not set by pair style in ewald/n");
double eps_i, sigma_i, sigma_n, *bi = B = new double[7*n+7];
double c[7] = {
1.0, sqrt(6.0), sqrt(15.0), sqrt(20.0), sqrt(15.0), sqrt(6.0), 1.0};
for (int i=0; i<=n; ++i) {
eps_i = sqrt(epsilon[i][i]);
sigma_i = sigma[i][i];
sigma_n = 1.0;
for (int j=0; j<7; ++j) {
*(bi++) = sigma_n*eps_i*c[j]; sigma_n *= sigma_i;
}
}
}
}
void EwaldN::init_coeff_sums() // sums based on atoms
{
if (sums) return; // calculated only once
sums = 1;
Sum sum_local[EWALD_MAX_NSUMS];
memset(sum_local, 0, EWALD_MAX_NSUMS*sizeof(Sum));
if (function[0]) { // 1/r
double *q = atom->q, *qn = q+atom->nlocal;
for (double *i=q; i<qn; ++i) {
sum_local[0].x += i[0]; sum_local[0].x2 += i[0]*i[0]; }
}
if (function[1]) { // geometric 1/r^6
int *type = atom->type, *ntype = type+atom->nlocal;
for (int *i=type; i<ntype; ++i) {
sum_local[1].x += B[i[0]]; sum_local[1].x2 += B[i[0]]*B[i[0]]; }
}
if (function[2]) { // aritmetic 1/r^6
double *bi;
int *type = atom->type, *ntype = type+atom->nlocal;
for (int *i=type; i<ntype; ++i) {
bi = B+7*i[0];
sum_local[2].x2 += bi[0]*bi[6];
for (int k=2; k<9; ++k) sum_local[k].x += *(bi++);
}
}
if (function[3]) { // dipole
int *type = atom->type, *ntype = type+atom->nlocal;
double *dipole = atom->dipole;
for (int *i=type; i<ntype; ++i)
sum_local[9].x2 += dipole[i[0]]*dipole[i[0]];
}
MPI_Allreduce(sum_local, sum, 2*EWALD_MAX_NSUMS, MPI_DOUBLE, MPI_SUM, world);
}
void EwaldN::init_self()
{
double g1 = g_ewald, g2 = g1*g1, g3 = g1*g2;
memset(energy_self, 0, EWALD_NFUNCS*sizeof(double)); // self energy
memset(virial_self, 0, EWALD_NFUNCS*sizeof(double));
if (function[0]) { // 1/r
virial_self[0] = -0.5*M_PI*qqrd2e/(g2*volume)*sum[0].x*sum[0].x;
energy_self[0] = sum[0].x2*qqrd2e*g1/sqrt(M_PI)-virial_self[0];
}
if (function[1]) { // geometric 1/r^6
virial_self[1] = M_PI*sqrt(M_PI)*g3/(6.0*volume)*sum[1].x*sum[1].x;
energy_self[1] = -sum[1].x2*g3*g3/12.0+virial_self[1];
}
if (function[2]) { // arithmetic 1/r^6
virial_self[2] = M_PI*sqrt(M_PI)*g3/(48.0*volume)*(sum[2].x*sum[8].x+
sum[3].x*sum[7].x+sum[4].x*sum[6].x+0.5*sum[5].x*sum[5].x);
energy_self[2] = -sum[2].x2*g3*g3/3.0+virial_self[2];
}
if (function[3]) { // dipole
virial_self[3] = 0; // in surface
energy_self[3] = sum[9].x2*mumurd2e*2.0*g3/3.0/sqrt(M_PI)-virial_self[3];
}
}
/* ----------------------------------------------------------------------
compute the EwaldN long-range force, energy, virial
------------------------------------------------------------------------- */
void EwaldN::compute(int eflag, int vflag)
{
if (!nbox) return;
reallocate_atoms();
compute_ek();
compute_force();
compute_surface();
compute_energy(eflag);
compute_virial(vflag);
}
void EwaldN::compute_ek()
{
cvector *ekr = ekr_local;
int lbytes = (2*nbox+1)*sizeof(cvector);
hvector *h;
kvector *k, *nk = kvec+nkvec;
cvector z1, z[2*nbox+1], *zx, *zy, *zz, *zn = z+2*nbox;
complex *cek, zxyz, zxy, cx;
vector mui;
double *x = atom->x[0], *xn = x+3*atom->nlocal, *q = atom->q, qi, bi, ci[7];
double *dipole = atom->dipole, *mu = atom->mu ? atom->mu[0] : NULL;
int i, kx, ky, n = nkvec*nsums, *type = atom->type, tri = domain->triclinic;
int func[EWALD_NFUNCS];
memcpy(func, function, EWALD_NFUNCS*sizeof(int));
memset(cek_local, 0, n*sizeof(complex)); // reset sums
while (x<xn) {
zx = (zy = (zz = z+nbox)+1)-2;
C_SET(zz->x, 1, 0); C_SET(zz->y, 1, 0); C_SET(zz->z, 1, 0); // z[0]
if (tri) { // triclinic z[1]
C_ANGLE(z1.x, unit[0]*x[0]+unit[5]*x[1]+unit[4]*x[2]);
C_ANGLE(z1.y, unit[1]*x[1]+unit[3]*x[2]);
C_ANGLE(z1.z, x[2]*unit[2]); x += 3;
}
else { // orthogonal z[1]
C_ANGLE(z1.x, *(x++)*unit[0]);
C_ANGLE(z1.y, *(x++)*unit[1]);
C_ANGLE(z1.z, *(x++)*unit[2]);
}
for (; zz<zn; --zx, ++zy, ++zz) { // set up z[k]=e^(ik.r)
C_RMULT(zy->x, zz->x, z1.x); // 3D k-vector
C_RMULT(zy->y, zz->y, z1.y); C_CONJ(zx->y, zy->y);
C_RMULT(zy->z, zz->z, z1.z); C_CONJ(zx->z, zy->z);
}
kx = ky = -1;
cek = cek_local;
if (func[0]) qi = *(q++);
if (func[1]) bi = B[*type];
if (func[2]) memcpy(ci, B+7*type[0], 7*sizeof(double));
if (func[3]) {
memcpy(mui, mu, sizeof(vector)); mu += 3;
vec_scalar_mult(mui, dipole[*type]);
h = hvec;
}
for (k=kvec; k<nk; ++k) { // compute rho(k)
if (ky!=k->y) { // based on order in
if (kx!=k->x) cx = z[kx = k->x].x; // reallocate
C_RMULT(zxy, z[ky = k->y].y, cx);
}
C_RMULT(zxyz, z[k->z].z, zxy);
if (func[0]) {
cek->re += zxyz.re*qi; (cek++)->im += zxyz.im*qi;
}
if (func[1]) {
cek->re += zxyz.re*bi; (cek++)->im += zxyz.im*bi;
}
if (func[2]) for (i=0; i<7; ++i) {
cek->re += zxyz.re*ci[i]; (cek++)->im += zxyz.im*ci[i];
}
if (func[3]) {
register double muk = mui[0]*h->x+mui[1]*h->y+mui[2]*h->z; ++h;
cek->re += zxyz.re*muk; (cek++)->im += zxyz.im*muk;
}
}
ekr = (cvector *) ((char *) memcpy(ekr, z, lbytes)+lbytes);
++type;
}
MPI_Allreduce(cek_local, cek_global, 2*n, MPI_DOUBLE, MPI_SUM, world);
}
void EwaldN::compute_force()
{
kvector *k;
hvector *h, *nh;
cvector *z = ekr_local;
vector sum[EWALD_MAX_NSUMS], mui;
complex *cek, zc, zx, zxy;
double *f = atom->f[0], *fn = f+3*atom->nlocal, *q = atom->q, *t = NULL;
double *dipole = atom->dipole, *mu = atom->mu ? atom->mu[0] : NULL;
double *ke, c[EWALD_NFUNCS] = {
8.0*M_PI*qqrd2e/volume, 2.0*M_PI*sqrt(M_PI)/(12.0*volume),
2.0*M_PI*sqrt(M_PI)/(192.0*volume), 8.0*M_PI*mumurd2e/volume};
double kt = 4.0*pow(g_ewald, 3.0)/3.0/sqrt(M_PI)/c[3];
int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type;
int func[EWALD_NFUNCS];
if (atom->torque) t = atom->torque[0];
memcpy(func, function, EWALD_NFUNCS*sizeof(int));
memset(sum, 0, EWALD_MAX_NSUMS*sizeof(vector)); // fj = -dE/dr =
for (; f<fn; f+=3) { // -i*qj*fac*
k = kvec; // Sum[conj(d)-d]
kx = ky = -1; // d = k*conj(ekj)*ek
ke = kenergy;
cek = cek_global;
memset(sum, 0, EWALD_MAX_NSUMS*sizeof(vector));
if (func[3]) {
register double di = dipole[*type]*c[3];
mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[1]; mui[2] = di*(mu++)[2];
}
for (nh = (h = hvec)+nkvec; h<nh; ++h, ++k) {
if (ky!=k->y) { // based on order in
if (kx!=k->x) zx = z[kx = k->x].x; // reallocate
C_RMULT(zxy, z[ky = k->y].y, zx);
}
C_CRMULT(zc, z[k->z].z, zxy);
if (func[0]) { // 1/r
register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek;
sum[0][0] += h->x*im; sum[0][1] += h->y*im; sum[0][2] += h->z*im;
}
if (func[1]) { // geometric 1/r^6
register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek;
sum[1][0] += h->x*im; sum[1][1] += h->y*im; sum[1][2] += h->z*im;
}
if (func[2]) { // arithmetic 1/r^6
register double im, c = *(ke++);
for (i=2; i<9; ++i) {
im = c*(zc.im*cek->re+cek->im*zc.re); ++cek;
sum[i][0] += h->x*im; sum[i][1] += h->y*im; sum[i][2] += h->z*im;
}
}
if (func[3]) { // dipole
register double im = *(ke++)*(zc.im*cek->re+
cek->im*zc.re)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z); ++cek;
sum[9][0] += h->x*im; sum[9][1] += h->y*im; sum[9][2] += h->z*im;
}
}
if (func[0]) { // 1/r
register double qi = *(q++)*c[0];
f[0] -= sum[0][0]*qi; f[1] -= sum[0][1]*qi; f[2] -= sum[0][2]*qi;
}
if (func[1]) { // geometric 1/r^6
register double bi = B[*type]*c[1];
f[0] -= sum[1][0]*bi; f[1] -= sum[1][1]*bi; f[2] -= sum[1][2]*bi;
}
if (func[2]) { // arithmetic 1/r^6
register double *bi = B+7*type[0]+7;
for (i=2; i<9; ++i) {
register double c2 = (--bi)[0]*c[2];
f[0] -= sum[i][0]*c2; f[1] -= sum[i][1]*c2; f[2] -= sum[i][2]*c2;
}
}
if (func[3]) { // dipole
f[0] -= sum[9][0]; f[1] -= sum[9][1]; f[2] -= sum[9][2];
*(t++) -= mui[1]*sum[0][2]+mui[2]*sum[0][1]-mui[0]*kt; // torque
*(t++) -= mui[2]*sum[0][0]+mui[0]*sum[0][2]-mui[1]*kt;
*(t++) -= mui[0]*sum[0][1]+mui[1]*sum[0][0]-mui[2]*kt;
}
z = (cvector *) ((char *) z+lbytes);
++type;
}
}
void EwaldN::compute_surface()
{
if (!function[3]) return;
vector sum_local = VECTOR_NULL, sum_total;
double *mu = atom->mu ? atom->mu[0] : NULL, *dipole = atom->dipole;
int *type = atom->type, *ntype = type+atom->nlocal;
for (int *i=type; i<ntype; ++i) {
register double di = dipole[i[0]];
sum_local[0] += di*(mu++)[0];
sum_local[1] += di*(mu++)[1];
sum_local[2] += di*(mu++)[2];
}
MPI_Allreduce(sum_local, sum_total, 3, MPI_DOUBLE, MPI_SUM, world);
energy_self[3] += virial_self[3];
virial_self[3] =
mumurd2e*(2.0*M_PI*vec_dot(sum_total,sum_total)/(2.0*dielectric+1)/volume);
energy_self[3] -= virial_self[3];
}
void EwaldN::compute_energy(int eflag)
{
energy = 0.0;
if (!eflag) return;
complex *cek = cek_global;
double *ke = kenergy;
double c[EWALD_NFUNCS] = {
4.0*M_PI*qqrd2e/volume, 2.0*M_PI*sqrt(M_PI)/(24.0*volume),
2.0*M_PI*sqrt(M_PI)/(192.0*volume), 4.0*M_PI*mumurd2e/volume};
double sum[EWALD_NFUNCS];
int func[EWALD_NFUNCS];
memcpy(func, function, EWALD_NFUNCS*sizeof(int));
memset(sum, 0, EWALD_NFUNCS*sizeof(double)); // reset sums
for (int k=0; k<nkvec; ++k) { // sum over k vectors
if (func[0]) { // 1/r
sum[0] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; }
if (func[1]) { // geometric 1/r^6
sum[1] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; }
if (func[2]) { // arithmetic 1/r^6
register double r =
(cek[0].re*cek[6].re+cek[0].im*cek[6].im)+
(cek[1].re*cek[5].re+cek[1].im*cek[5].im)+
(cek[2].re*cek[4].re+cek[2].im*cek[4].im)+
0.5*(cek[3].re*cek[3].re+cek[3].im*cek[3].im); cek += 7;
sum[2] += *(ke++)*r;
}
if (func[3]) { // dipole
sum[3] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; }
}
for (int k=0; k<EWALD_NFUNCS; ++k) energy += c[k]*sum[k]-energy_self[k];
if (slabflag) compute_slabcorr(eflag);
}
#define swap(a, b) { register double t = a; a= b; b = t; }
void EwaldN::compute_virial(int vflag)
{
memset(virial, 0, sizeof(shape));
if (!vflag) return;
complex *cek = cek_global;
double *kv = kvirial;
double c[EWALD_NFUNCS] = {
4.0*M_PI*qqrd2e/volume, 2.0*M_PI*sqrt(M_PI)/(24.0*volume),
2.0*M_PI*sqrt(M_PI)/(192.0*volume), 4.0*M_PI*mumurd2e/volume};
shape sum[EWALD_NFUNCS];
int func[EWALD_NFUNCS];
memcpy(func, function, EWALD_NFUNCS*sizeof(int));
memset(sum, 0, EWALD_NFUNCS*sizeof(shape));
for (int k=0; k<nkvec; ++k) { // sum over k vectors
if (func[0]) { // 1/r
register double r = cek->re*cek->re+cek->im*cek->im; ++cek;
sum[0][0] += *(kv++)*r; sum[0][1] += *(kv++)*r; sum[0][2] += *(kv++)*r;
sum[0][3] += *(kv++)*r; sum[0][4] += *(kv++)*r; sum[0][5] += *(kv++)*r;
}
if (func[1]) { // geometric 1/r^6
register double r = cek->re*cek->re+cek->im*cek->im; ++cek;
sum[1][0] += *(kv++)*r; sum[1][1] += *(kv++)*r; sum[1][2] += *(kv++)*r;
sum[1][3] += *(kv++)*r; sum[1][4] += *(kv++)*r; sum[1][5] += *(kv++)*r;
}
if (func[2]) { // arithmetic 1/r^6
register double r =
(cek[0].re*cek[6].re+cek[0].im*cek[6].im)+
(cek[1].re*cek[5].re+cek[1].im*cek[5].im)+
(cek[2].re*cek[4].re+cek[2].im*cek[4].im)+
0.5*(cek[3].re*cek[3].re+cek[3].im*cek[3].im); cek += 7;
sum[2][0] += *(kv++)*r; sum[2][1] += *(kv++)*r; sum[2][2] += *(kv++)*r;
sum[2][3] += *(kv++)*r; sum[2][4] += *(kv++)*r; sum[2][5] += *(kv++)*r;
}
if (func[3]) {
register double r = cek->re*cek->re+cek->im*cek->im; ++cek;
sum[3][0] += *(kv++)*r; sum[3][1] += *(kv++)*r; sum[3][2] += *(kv++)*r;
sum[3][3] += *(kv++)*r; sum[3][4] += *(kv++)*r; sum[3][5] += *(kv++)*r;
}
}
for (int k=0; k<EWALD_NFUNCS; ++k)
if (func[k]) {
shape self = {virial_self[k], virial_self[k], virial_self[k], 0, 0, 0};
shape_scalar_mult(sum[k], c[k]);
shape_add(virial, sum[k]);
shape_subtr(virial, self);
}
}
/* ----------------------------------------------------------------------
Slab-geometry correction term to dampen inter-slab interactions between
periodically repeating slabs. Yields good approximation to 2-D EwaldN if
adequate empty space is left between repeating slabs (J. Chem. Phys.
111, 3155). Slabs defined here to be parallel to the xy plane.
------------------------------------------------------------------------- */
void EwaldN::compute_slabcorr(int eflag)
{
// compute local contribution to global dipole moment
double *q = atom->q;
double *x = atom->x[0]-1, *xn = x+3*atom->nlocal-1;
double dipole = 0.0, dipole_all;
while ((x+=3)<xn) dipole += *x * *(q++);
MPI_Allreduce(&dipole, &dipole_all, 1, MPI_DOUBLE, MPI_SUM, world);
double ffact = -4.0*M_PI*qqrd2e*dipole_all/volume; // force correction
double *f = atom->f[0]-1, *fn = f+3*atom->nlocal-3;
q = atom->q;
while ((f+=3)<fn) *f += ffact* *(q++);
if (eflag) // energy correction
energy += qqrd2e*(2.0*M_PI*dipole_all*dipole_all/volume);
}

75
src/USER-EWALDN/ewald_n.h Normal file
View File

@ -0,0 +1,75 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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.
------------------------------------------------------------------------- */
#ifndef EWALD_N_H
#define EWALD_N_H
#include "kspace.h"
#include "math_complex.h"
#define EWALD_NORDER 6
#define EWALD_NFUNCS 4
#define EWALD_MAX_NSUMS 10
#define EWALD_NSUMS {1, 1, 7, 1}
namespace LAMMPS_NS {
typedef struct cvector { complex x, y, z; } cvector;
typedef struct hvector { double x, y, z; } hvector;
typedef struct kvector { long x, y, z; } kvector;
class EwaldN : public KSpace {
public:
EwaldN(class LAMMPS *, int, char **);
~EwaldN();
void init();
void setup();
void compute(int, int);
int memory_usage() { return bytes; }
private:
double unit[6];
int function[EWALD_NFUNCS], first_output;
int nkvec, nkvec_max, nevec, nevec_max,
nbox, nfunctions, nsums, bytes, sums;
double precision, g2_max;
double *kenergy, energy_self[EWALD_NFUNCS];
double *kvirial, virial_self[EWALD_NFUNCS];
cvector *ekr_local;
hvector *hvec;
kvector *kvec;
double qqrd2e, mumurd2e, dielectric, *B, volume;
struct Sum { double x, x2; } sum[EWALD_MAX_NSUMS];
complex *cek_local, *cek_global;
void reallocate();
void reallocate_atoms();
void deallocate();
void coefficients();
void init_coeffs();
void init_coeff_sums();
void init_self();
void compute_ek();
void compute_force();
void compute_surface();
void compute_energy(int);
void compute_virial(int);
void compute_slabcorr(int);
};
}
#endif

73
src/USER-EWALDN/math_complex.h Normal file
View File

@ -0,0 +1,73 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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: Pieter J. in 't Veld (SNL)
------------------------------------------------------------------------- */
#ifndef MATH_COMPLEX_H
#define MATH_COMPLEX_H
#define COMPLEX_NULL {0, 0}
namespace LAMMPS_NS {
typedef struct complex {
double re, im; } complex;
}
#define C_MULT(d, x, y) { \
d.re = x.re*y.re-x.im*y.im; \
d.im = x.re*y.im+x.im*y.re; }
#define C_RMULT(d, x, y) { \
register complex t = x; \
d.re = t.re*y.re-t.im*y.im; \
d.im = t.re*y.im+t.im*y.re; }
#define C_CRMULT(d, x, y) { \
register complex t = x; \
d.re = t.re*y.re-t.im*y.im; \
d.im = -t.re*y.im-t.im*y.re; }
#define C_SMULT(d, x, y) { \
d.re = x.re*y; \
d.im = x.im*y; }
#define C_ADD(d, x, y) { \
d.re = x.re+y.re; \
d.im = x.im+y.im; }
#define C_SUBTR(d, x, y) { \
d.re = x.re-y.re; \
d.im = x.im-y.im; }
#define C_CONJ(d, x) { \
d.re = x.re; \
d.im = -x.im; }
#define C_SET(d, x, y) { \
d.re = x; \
d.im = y; }
#define C_ANGLE(d, angle) { \
register double a = angle; \
d.re = cos(a); \
d.im = sin(a); }
#define C_COPY(d, x) { \
memcpy(&d, &x, sizeof(complex)); }
#endif

447
src/USER-EWALDN/math_vector.h Normal file
View File

@ -0,0 +1,447 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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: Pieter J. in 't Veld (SNL)
------------------------------------------------------------------------- */
#ifndef MATH_VECTOR_H
#define MATH_VECTOR_H
#include "math.h"
#define VECTOR_NULL {0, 0, 0}
#define SHAPE_NULL {0, 0, 0, 0, 0, 0}
#define FORM_NULL {0, 0, 0, 0, 0, 0}
#define MATRIX_NULL {VECTOR_NULL, VECTOR_NULL, VECTOR_NULL}
#define VECTOR4_NULL {0, 0, 0, 0}
#define QUATERNION_NULL {0, 0, 0, 0}
#define FORM4_NULL {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
#define FZERO 1e-15
#define fzero(x) (((x)>-FZERO) && ((x)<FZERO))
namespace LAMMPS_NS {
typedef double vector[3]; // 0:x 1:y 2:z
typedef int lvector[3];
typedef double shape[6]; // 0:xx 1:yy 2:zz 3:zy 4:zx 5:yx
typedef int lshape[6]; // xy=0 xz=0 yz=0;
typedef double form[6]; // 0:xx 1:yy 2:zz 3:zy 4:zx 5:yx
typedef int lform[6]; // xy=yx xz=zx yz=zy;
typedef vector matrix[3]; // 00:xx 11:yy 22:zz 21:zy 20:zx 10:yx
typedef lvector lmatrix[3];
typedef double vector4[4]; // 4D vector
typedef double form4[10]; // 0:00 1:11 2:22 3:33 4:32
// 5:31 6:30 7:21 8:20 9:10
// 01=10 02=20 03=30 etc
typedef double quaternion[4]; // quaternion
// vector operators
inline double vec_dot(vector &a, vector &b) { // a.b
return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; }
inline void vec_null(vector &dest) {
memset(dest, 0, sizeof(vector)); }
inline void vec_neg(vector &dest) { // -a
dest[0] = -dest[0];
dest[1] = -dest[1];
dest[2] = -dest[2]; }
inline void vec_norm(vector &dest) { // a/|a|
register double f = sqrt(vec_dot(dest, dest));
dest[0] /= f;
dest[1] /= f;
dest[2] /= f; }
inline void vec_add(vector &dest, vector &src) { // a+b
dest[0] += src[0];
dest[1] += src[1];
dest[2] += src[2]; }
inline void vec_subtr(vector &dest, vector &src) { // a-b
dest[0] -= src[0];
dest[1] -= src[1];
dest[2] -= src[2]; }
inline void vec_mult(vector &dest, vector &src) { // a*b
dest[0] *= src[0];
dest[1] *= src[1];
dest[2] *= src[2]; }
inline void vec_div(vector &dest, vector &src) { // a/b
dest[0] /= src[0];
dest[1] /= src[1];
dest[2] /= src[2]; }
inline void vec_cross(vector &dest, vector &src) { // a x b
vector t = {
dest[1]*src[2]-dest[2]*src[1],
dest[2]*src[0]-dest[0]*src[2],
dest[0]*src[1]-dest[1]*src[0]};
memcpy(dest, t, sizeof(vector)); }
inline void vec_scalar_mult(vector &dest, double f) { // f*a
dest[0] *= f;
dest[1] *= f;
dest[2] *= f; }
inline void vec_to_lvec(lvector &dest, vector &src) {
dest[0] = (int) src[0];
dest[1] = (int) src[1];
dest[2] = (int) src[2]; }
inline void lvec_to_vec(vector &dest, lvector &src) {
dest[0] = (double) src[0];
dest[1] = (double) src[1];
dest[2] = (double) src[2]; }
// shape operators
inline double shape_det(shape &s) {
return s[0]*s[1]*s[2]; }
inline void shape_null(shape &dest) {
memset(dest, 0, sizeof(shape)); }
inline void shape_unit(shape &dest) {
memset(dest, 0, sizeof(shape));
dest[0] = dest[1] = dest[2] = 1.0; }
inline void shape_add(shape &dest, shape &src) { // h_a+h_b
dest[0] += src[0]; dest[1] += src[1]; dest[2] += src[2];
dest[3] += src[3]; dest[4] += src[4]; dest[5] += src[5]; }
inline void shape_subtr(shape &dest, shape &src) { // h_a-h_b
dest[0] -= src[0]; dest[1] -= src[1]; dest[2] -= src[2];
dest[3] -= src[3]; dest[4] -= src[4]; dest[5] -= src[5]; }
inline void shape_inv(shape &h_inv, shape &h) { // h^-1
h_inv[0] = 1.0/h[0]; h_inv[1] = 1.0/h[1]; h_inv[2] = 1.0/h[2];
h_inv[3] = -h[3]/(h[1]*h[2]);
h_inv[4] = (h[3]*h[5]-h[1]*h[4])/(h[0]*h[1]*h[2]);
h_inv[5] = -h[5]/(h[0]*h[1]); }
inline void shape_dot(shape &dest, shape &src) { // h_a.h_b
dest[3] = dest[1]*src[3]+dest[3]*src[2];
dest[4] = dest[0]*src[4]+dest[5]*src[3]+dest[4]*src[2];
dest[5] = dest[0]*src[5]+dest[5]*src[1];
dest[0] *= src[0]; dest[1] *= src[1]; dest[2] *= src[2]; }
inline void shape_scalar_mult(shape &dest, double f) { // f*h
dest[0] *= f; dest[1] *= f; dest[2] *= f;
dest[3] *= f; dest[4] *= f; dest[5] *= f; }
inline void shape_vec_dot(vector &dest, vector &src, shape &h) {// h.a
dest[0] = h[0]*src[0]+h[5]*src[1]+h[4]*src[2];
dest[1] = h[1]*src[1]+h[3]*src[2];
dest[2] = h[2]*src[2]; }
inline void vec_shape_dot(vector &dest, shape &h, vector &src) {// a.h
dest[2] = h[4]*src[0]+h[3]*src[1]+h[2]*src[2];
dest[1] = h[5]*src[0]+h[1]*src[1];
dest[0] = h[0]*src[0]; }
inline void shape_to_matrix(matrix &dest, shape &h) { // m = h
dest[0][0] = h[0]; dest[1][0] = h[5]; dest[2][0] = h[4];
dest[0][1] = 0.0; dest[1][1] = h[1]; dest[2][1] = h[3];
dest[0][2] = 0.0; dest[1][2] = 0.0; dest[2][2] = h[2]; }
inline void shape_to_lshape(lshape &dest, shape &src) {
dest[0] = (long)src[0]; dest[1] = (long)src[1]; dest[2] = (long)src[2];
dest[3] = (long)src[3]; dest[4] = (long)src[4]; dest[5] = (long)src[5]; }
inline void lshape_to_shape(shape &dest, lshape &src) {
dest[0] = (double)src[0]; dest[1] = (double)src[1]; dest[2] = (double)src[2];
dest[3] = (double)src[3]; dest[4] = (double)src[4]; dest[5] = (double)src[5];}
// form operators
inline double form_det(form &m) { // |m|
return m[0]*(m[1]*m[2]-m[3]*m[3])+
m[4]*(2.0*m[3]*m[5]-m[1]*m[4])-m[2]*m[5]*m[5]; }
inline void form_null(form &dest) {
memset(dest, 0, sizeof(form)); }
inline void form_unit(form &dest) {
memset(dest, 0, sizeof(form));
dest[0] = dest[1] = dest[2] = 1.0; }
inline void form_add(form &dest, form &src) { // m_a+m_b
dest[0] += src[0]; dest[1] += src[1]; dest[2] += src[2];
dest[3] += src[3]; dest[4] += src[4]; dest[5] += src[5]; }
inline void form_subtr(shape &dest, form &src) { // m_a-m_b
dest[0] -= src[0]; dest[1] -= src[1]; dest[2] -= src[2];
dest[3] -= src[3]; dest[4] -= src[4]; dest[5] -= src[5]; }
inline int form_inv(form &m_inv, form &m) { // m^-1
register double det = form_det(m);
if (fzero(det)) return 0;
m_inv[0] = (m[1]*m[2]-m[3]*m[3])/det;
m_inv[1] = (m[0]*m[2]-m[4]*m[4])/det;
m_inv[2] = (m[0]*m[1]-m[5]*m[5])/det;
m_inv[3] = (m[4]*m[5]-m[0]*m[3])/det;
m_inv[4] = (m[3]*m[5]-m[1]*m[4])/det;
m_inv[5] = (m[3]*m[4]-m[2]*m[5])/det;
return 1; }
inline void form_dot(form &dest, form &src) { // m_a.m_b
form m;
memcpy(m, dest, sizeof(form));
dest[0] = m[0]*src[0]+m[4]*src[4]+m[5]*src[5];
dest[1] = m[1]*src[1]+m[3]*src[3]+m[5]*src[5];
dest[2] = m[2]*src[2]+m[3]*src[3]+m[4]*src[4];
dest[3] = m[3]*src[2]+m[1]*src[3]+m[5]*src[4];
dest[4] = m[4]*src[2]+m[5]*src[3]+m[0]*src[4];
dest[5] = m[5]*src[1]+m[4]*src[3]+m[0]*src[5]; }
inline void form_vec_dot(vector &dest, form &m) { // m.a
vector a;
memcpy(a, dest, sizeof(vector));
dest[0] = m[0]*a[0]+m[5]*a[1]+m[4]*a[2];
dest[1] = m[5]*a[0]+m[1]*a[1]+m[3]*a[2];
dest[2] = m[4]*a[0]+m[3]*a[1]+m[2]*a[2]; }
inline void form_to_lform(lform &dest, form &src) {
dest[0] = (long)src[0]; dest[1] = (long)src[1]; dest[2] = (long)src[2];
dest[3] = (long)src[3]; dest[4] = (long)src[4]; dest[5] = (long)src[5]; }
inline void lform_to_form(form &dest, lform &src) {
dest[0] = (double)src[0]; dest[1] = (double)src[1]; dest[2] = (double)src[2];
dest[3] = (double)src[3]; dest[4] = (double)src[4]; dest[5] = (double)src[5];}
// matrix operators
inline double matrix_det(matrix &m) { // |m|
vector axb;
memcpy(&axb, m[0], sizeof(vector));
vec_cross(axb, m[1]);
return vec_dot(axb, m[2]); }
inline void matrix_null(matrix &dest) {
memset(dest, 0, sizeof(dest)); }
inline void matrix_unit(matrix &dest) {
memset(dest, 0, sizeof(dest));
dest[0][0] = dest[1][1] = dest[2][2] = 1.0; }
inline void matrix_scalar_mult(matrix &dest, double f) { // f*m
vec_scalar_mult(dest[0], f);
vec_scalar_mult(dest[1], f);
vec_scalar_mult(dest[2], f); }
inline void matrix_trans(matrix &dest) { // m^t
double f = dest[0][1]; dest[0][1] = dest[1][0]; dest[1][0] = f;
f = dest[0][2]; dest[0][2] = dest[2][0]; dest[2][0] = f;
f = dest[1][2]; dest[1][2] = dest[2][1]; dest[2][1] = f; }
inline int matrix_inv(matrix &dest, matrix &src) { // m^-1
double f = matrix_det(src);
if (fzero(f)) return 0; // singular matrix
memcpy(dest[0], src[1], sizeof(vector));
memcpy(dest[1], src[2], sizeof(vector));
memcpy(dest[2], src[0], sizeof(vector));
vec_cross(dest[0], src[2]);
vec_cross(dest[1], src[0]);
vec_cross(dest[2], src[1]);
matrix_scalar_mult(dest, 1.0/f);
matrix_trans(dest);
return 0; }
inline void matrix_vec_dot(vector &dest, vector &src, matrix &m) { // m.a
dest[0] = m[0][0]*src[0]+m[1][0]*src[1]+m[2][0]*src[2];
dest[1] = m[0][1]*src[0]+m[1][1]*src[1]+m[2][1]*src[2];
dest[2] = m[0][2]*src[0]+m[1][2]*src[1]+m[2][2]*src[2]; }
inline void matrix_to_shape(shape &dest, matrix &src) { // h = m
dest[0] = src[0][0]; dest[1] = src[1][1]; dest[2] = src[2][2];
dest[3] = src[2][1]; dest[4] = src[2][0]; dest[5] = src[1][0]; }
inline void matrix_to_lmatrix(lmatrix &dest, matrix &src) {
vec_to_lvec(dest[0], src[0]);
vec_to_lvec(dest[1], src[1]);
vec_to_lvec(dest[2], src[2]); }
inline void lmatrix_to_matrix(matrix &dest, lmatrix &src) {
lvec_to_vec(dest[0], src[0]);
lvec_to_vec(dest[1], src[1]);
lvec_to_vec(dest[2], src[2]); }
// quaternion operators
inline double quat_dot(quaternion &p, quaternion &q) { // p.q
return p[0]*q[0]+p[1]*q[1]+p[2]*q[2]+p[3]*q[3];
}
inline void quat_norm(quaternion &q) { // q = q/|q|
double f = sqrt(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]+q[3]*q[3]);
q[0] /= f; q[1] /= f; q[2] /= f; q[3] /= f;
}
inline void quat_conj(quaternion &q) { // q = conj(q)
q[1] = -q[1]; q[2] = -q[2]; q[3] = -q[3];
}
inline void quat_mult(quaternion &dest, quaternion &src) { // dest *= src
quaternion q;
memcpy(q, dest, sizeof(quaternion));
dest[0] = src[0]*q[3]+src[1]*q[2]-src[2]*q[1]+src[3]*q[0];
dest[1] = -src[0]*q[2]+src[1]*q[3]+src[2]*q[0]+src[3]*q[1];
dest[2] = src[0]*q[1]-src[1]*q[0]+src[2]*q[3]+src[3]*q[2];
dest[3] = -src[0]*q[0]-src[1]*q[1]-src[2]*q[2]+src[3]*q[3];
}
inline void quat_div(quaternion &dest, quaternion &src) { // dest /= src
quaternion q;
memcpy(q, dest, sizeof(quaternion));
dest[0] = src[0]*q[3]-src[1]*q[2]+src[2]*q[1]-src[3]*q[0];
dest[1] = -src[0]*q[2]-src[1]*q[3]-src[2]*q[0]-src[3]*q[1];
dest[2] = src[0]*q[1]+src[1]*q[0]-src[2]*q[3]-src[3]*q[2];
dest[3] = -src[0]*q[0]+src[1]*q[1]+src[2]*q[2]-src[3]*q[3];
}
// dest = q*src*conj(q)
inline void quat_vec_rot(vector &dest, vector &src, quaternion &q) {
quaternion aa={q[0]*q[0], q[1]*q[1], q[2]*q[2], q[3]*q[3]};
form ab={q[0]*q[1], q[0]*q[2], q[0]*q[3], q[1]*q[2], q[1]*q[3], q[2]*q[3]};
dest[0] = (aa[0]+aa[1]-aa[2]-aa[3])*src[0]+
((ab[3]-ab[2])*src[1]+(ab[1]+ab[4])*src[2])*2.0;
dest[1] = (aa[0]-aa[1]+aa[2]-aa[3])*src[1]+
((ab[2]+ab[3])*src[0]+(ab[5]-ab[0])*src[2])*2.0;
dest[2] = (aa[0]-aa[1]-aa[2]+aa[3])*src[2]+
((ab[4]-ab[1])*src[0]+(ab[0]+ab[5])*src[1])*2.0;
}
// vector4 operators
inline void vec4_null(vector4 &dest) {
memset(dest, 0, sizeof(vector4));
}
inline double vec4_dot(vector4 &a, vector4 &b) {
return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]+a[3]*b[3]; }
// form operators
inline void form4_null(form4 &dest) {
memset(dest, 0, sizeof(form4)); }
inline void form4_unit(form4 &dest) {
memset(dest, 0, sizeof(form4));
dest[0] = dest[1] = dest[2] = dest[3] = 1.0; }
inline double form4_det(form4 &m) {
register double f = m[6]*m[7]-m[5]*m[8];
return m[0]*(
m[1]*(m[2]*m[3]-m[4]*m[4])+
m[5]*(2.0*m[4]*m[7]-m[2]*m[5])-m[3]*m[7]*m[7])+f*f+
-m[1]*(m[2]*m[6]*m[6]+m[8]*(m[3]*m[8]-2.0*m[4]*m[6]))+
m[9]*(
2.0*(m[2]*m[5]*m[6]+m[3]*m[7]*m[8]-m[4]*(m[6]*m[7]+m[5]*m[8]))+
m[9]*(m[4]*m[4]-m[2]*m[3])); }
inline int form4_inv(form4 &m_inv, form4 &m) {
register double det = form4_det(m);
if (fzero(det)) return 0;
m_inv[0] = (m[1]*(m[2]*m[3]-m[4]*m[4])+
m[5]*(2.0*m[4]*m[7]-m[2]*m[5])-m[3]*m[7]*m[7])/det;
m_inv[1] = (m[0]*(m[2]*m[3]-m[4]*m[4])+
m[6]*(2.0*m[4]*m[8]-m[2]*m[6])-m[3]*m[8]*m[8])/det;
m_inv[2] = (m[0]*(m[1]*m[3]-m[5]*m[5])+
m[6]*(2.0*m[5]*m[9]-m[1]*m[6])-m[3]*m[9]*m[9])/det;
m_inv[3] = (m[0]*(m[1]*m[2]-m[7]*m[7])+
m[8]*(2.0*m[7]*m[9]-m[1]*m[8])-m[2]*m[9]*m[9])/det;
m_inv[4] = (m[0]*(m[5]*m[7]-m[1]*m[4])+m[1]*m[6]*m[8]+
m[9]*(m[4]*m[9]-m[6]*m[7]-m[5]*m[8]))/det;
m_inv[5] = (m[0]*(m[4]*m[7]-m[2]*m[5])+m[2]*m[6]*m[9]+
m[8]*(m[5]*m[8]-m[6]*m[7]-m[4]*m[9]))/det;
m_inv[6] = (m[1]*(m[4]*m[8]-m[2]*m[6])+m[2]*m[5]*m[9]+
m[7]*(m[6]*m[7]-m[5]*m[8]-m[4]*m[9]))/det;
m_inv[7] = (m[0]*(m[4]*m[5]-m[3]*m[7])+m[3]*m[8]*m[9]+
m[6]*(m[6]*m[7]-m[5]*m[8]-m[4]*m[9]))/det;
m_inv[8] = (m[1]*(m[4]*m[6]-m[3]*m[8])+m[3]*m[7]*m[9]+
m[5]*(m[5]*m[8]-m[6]*m[7]-m[4]*m[9]))/det;
m_inv[9] = (m[2]*(m[5]*m[6]-m[3]*m[9])+m[3]*m[7]*m[8]+
m[4]*(m[4]*m[9]-m[6]*m[7]-m[5]*m[8]))/det;
return 1; }
inline void form4_vec4_dot(vector4 &dest, form4 &m) {
vector4 a;
memcpy(a, dest, sizeof(vector4));
dest[0] = m[0]*a[0]+m[9]*a[1]+m[7]*a[2]+m[6]*a[3];
dest[1] = m[9]*a[0]+m[1]*a[1]+m[6]*a[2]+m[5]*a[3];
dest[2] = m[8]*a[0]+m[7]*a[1]+m[2]*a[2]+m[4]*a[3];
dest[3] = m[6]*a[0]+m[5]*a[1]+m[4]*a[2]+m[3]*a[3]; }
// square regression: y = eqn[0] + eqn[1]*x + eqn[2]*x*x
inline int regress2(vector &eqn, int order, double *x, double *y, int n) {
form xtx = FORM_NULL, xtx_inv;
vector xty = VECTOR_NULL;
double xn, xi, yi;
int i;
vec_null(eqn);
xtx[0] = n;
if ((order = order%2)<0) order = -order; // max: quad regress
if (order<1) xtx[1] = 1.0;
if (order<2) xtx[2] = 1.0;
for (i=0; i<n; ++i) {
xty[0] += (yi = y[i]);
if (order<1) continue;
xty[1] += yi*(xi = xn = x[i]); xtx[5] += xn; xtx[1] += (xn *= xi);
if (order<2) continue;
xty[2] += yi*xn; xtx[3] += (xn *= xi); xtx[2] += xn*xi;
}
if (order>1) xtx[4] = xtx[2];
if (!form_inv(xtx_inv, xtx)) return 0;
memcpy(eqn, xty, sizeof(vector));
form_vec_dot(eqn, xtx_inv);
return 1; }
// cubic regression: y = eqn[0] + eqn[1]*x + eqn[2]*x*x + eqn[3]*x*x*x
inline int regress3(vector4 &eqn, int order, double *x, double *y, int n) {
form4 xtx = FORM4_NULL, xtx_inv;
vector4 xty = VECTOR4_NULL;
double xn, xi, yi;
int i;
vec4_null(eqn);
xtx[0] = n;
if ((order = order%3)<0) order = -order; // max: cubic regress
if (order<1) xtx[1] = 1.0;
if (order<2) xtx[2] = 1.0;
if (order<3) xtx[3] = 1.0;
for (i=0; i<n; ++i) {
xty[0] += (yi = y[i]);
if (order<1) continue;
xty[1] += yi*(xi = xn = x[i]); xtx[9] += xn; xtx[1] += (xn *= xi);
if (order<2) continue;
xty[2] += yi*xn; xtx[7] += (xn *= xi); xtx[2] += xn*xi;
if (order<3) continue;
xty[3] += yi*xn; xtx[4] += (xn *= xi*xi); xtx[3] += xn*xi;
}
if (order>1) xtx[8] = xtx[1];
if (order>2) { xtx[6] = xtx[7]; xtx[5] = xtx[2]; }
if (!form4_inv(xtx_inv, xtx)) return 0;
memcpy(eqn, xty, sizeof(vector4));
form4_vec4_dot(eqn, xtx_inv);
return 1; }
}
#endif

1176
src/USER-EWALDN/pair_buck_coul.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

69
src/USER-EWALDN/pair_buck_coul.h Normal file
View File

@ -0,0 +1,69 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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.
------------------------------------------------------------------------- */
#ifndef PAIR_BUCK_COUL_H
#define PAIR_BUCK_COUL_H
#include "pair.h"
namespace LAMMPS_NS {
class PairBuckCoul : public Pair {
public:
double cut_coul;
PairBuckCoul(class LAMMPS *);
~PairBuckCoul();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
virtual void init_style();
void write_restart(FILE *);
void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual void single(int, int, int, int, double, double, double, int, One &);
virtual void *extract_ptr(char *);
virtual void extract_long(double *);
void compute_inner();
void compute_middle();
void compute_outer(int, int);
protected:
double cut_buck_global;
double **cut_buck, **cut_buck_read, **cut_bucksq;
double cut_coulsq;
double **buck_a_read, **buck_a, **buck_c_read, **buck_c;
double **buck1, **buck2, **buck_rho_read, **buck_rho, **rhoinv, **offset;
double *cut_respa;
double g_ewald;
int ewald_order, ewald_off;
double tabinnersq;
double *rtable, *drtable, *ftable, *dftable, *ctable, *dctable;
double *etable, *detable, *ptable, *dptable, *vtable, *dvtable;
int ncoulshiftbits, ncoulmask;
void options(char **arg, int order);
void allocate();
void init_tables();
void free_tables();
};
}
#endif

1159
src/USER-EWALDN/pair_lj_coul.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

68
src/USER-EWALDN/pair_lj_coul.h Normal file
View File

@ -0,0 +1,68 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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.
------------------------------------------------------------------------- */
#ifndef PAIR_LJ_COUL_H
#define PAIR_LJ_COUL_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJCoul : public Pair {
public:
double cut_coul;
PairLJCoul(class LAMMPS *);
virtual ~PairLJCoul();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
virtual void init_style();
void write_restart(FILE *);
void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual void single(int, int, int, int, double, double, double, int, One &);
virtual void *extract_ptr(char *);
virtual void extract_long(double *);
void compute_inner();
void compute_middle();
void compute_outer(int, int);
protected:
double cut_lj_global;
double **cut_lj, **cut_lj_read, **cut_ljsq;
double cut_coulsq;
double **epsilon_read, **epsilon, **sigma_read, **sigma;
double **lj1, **lj2, **lj3, **lj4, **offset;
double *cut_respa;
double g_ewald;
int ewald_order, ewald_off;
double tabinnersq;
double *rtable, *drtable, *ftable, *dftable, *ctable, *dctable;
double *etable, *detable, *ptable, *dptable, *vtable, *dvtable;
int ncoulshiftbits, ncoulmask;
void options(char **arg, int order);
void allocate();
void init_tables();
void free_tables();
};
}
#endif

30
src/USER-EWALDN/style_user_ewaldn.h Normal file
View File

@ -0,0 +1,30 @@
/* ----------------------------------------------------------------------
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 KSpaceInclude
#include "ewald_n.h"
#endif
#ifdef KSpaceClass
KSpaceStyle(ewald/n,EwaldN)
#endif
#ifdef PairInclude
#include "pair_buck_coul.h"
#include "pair_lj_coul.h"
#endif
#ifdef PairClass
PairStyle(buck/coul,PairBuckCoul)
PairStyle(lj/coul,PairLJCoul)
#endif

0
src/style_user_ackland.h Normal file
View File

0
src/style_user_ewaldn.h Normal file
View File

18
src/style_user_packages.h Normal file
View File

@ -0,0 +1,18 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
// style flies for user-contributed packages
// see the README files in individual user-package directories for details
#include "style_user_ackland.h"
#include "style_user_ewaldn.h"