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

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

This commit is contained in:
sjplimp 2014-09-05 23:16:15 +00:00
parent a6e4212209
commit af21769abe
12 changed files with 5202 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
src/Makefile
View File

@ -15,7 +15,7 @@ OBJ = $(SRC:.cpp=.o)
PACKAGE = asphere body class2 colloid dipole fld gpu granular kim \
kokkos kspace manybody mc meam misc molecule mpiio opt peri poems \
reax replica rigid shock srd voronoi xtc
reax replica rigid shock snap srd voronoi xtc
PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
user-cuda user-eff user-fep user-intel user-lb user-misc \
@ -38,8 +38,6 @@ help:
@echo ''
@echo 'make clean-all delete all object files'
@echo 'make clean-machine delete object files for one machine'
@echo 'make purge purge obsolete copies of package sources'
@echo 'make tar create lmp_src.tar.gz of src dir and packages'
@echo 'make makelib create Makefile.lib for static library build'
@echo 'make makeshlib create Makefile.shlib for shared library build'
@echo 'make makelist create Makefile.list used by old makes'
@ -48,6 +46,7 @@ help:
@echo 'make -f Makefile.list machine build LAMMPS from explicit list of files'
@echo 'make stubs build dummy MPI library in STUBS'
@echo 'make install-python install LAMMPS wrapper in Python'
@echo 'make tar create lmp_src.tar.gz of src dir and packages'
@echo ''
@echo 'make package list available packages'
@echo 'make package-status (ps) status of all packages'
@ -64,20 +63,29 @@ help:
@echo 'make package-update (pu) replace src files with updated package files'
@echo 'make package-overwrite replace package files with src files'
@echo 'make package-diff (pd) diff src files against package files'
@echo 'make package-purge purge obsolete copies of package sources'
@echo ''
@echo 'make machine build LAMMPS where machine is one of:'
@echo ''
@files="`ls MAKE/Makefile.*`"; \
for file in $$files; do head -1 $$file; done
@echo ''
@echo 'or one of:'
@echo ''
@files="`ls MAKE/MORE/Makefile.*`"; \
for file in $$files; do head -1 $$file; done
@echo ''
# Build the code
.DEFAULT:
@test -f MAKE/Makefile.$@
@test -f MAKE/Makefile.$@ -o -f MAKE/MORE/Makefile.$@
@if [ ! -d Obj_$@ ]; then mkdir Obj_$@; fi
@$(SHELL) Make.sh style
@cp MAKE/Makefile.$@ Obj_$@/Makefile
@if [ -f MAKE/Makefile.$@ ]; \
then cp MAKE/Makefile.$@ Obj_$@/Makefile; fi
@if [ -f MAKE/MORE/Makefile.$@ ]; \
then cp MAKE/MORE/Makefile.$@ Obj_$@/Makefile; fi
@if [ ! -e Makefile.package ]; \
then cp Makefile.package.empty Makefile.package; fi
@if [ ! -e Makefile.package.settings ]; \
@ -99,28 +107,6 @@ clean-all:
clean-%:
rm -rf Obj_$(@:clean-%=%)
purge: Purge.list
@echo 'Purging obsolete and auto-generated source files'
@for f in `grep -v '#' Purge.list` ; \
do test -f $$f && rm $$f && echo $$f || : ; \
done
# Create a tarball of src dir and packages
tar:
@cd STUBS; $(MAKE) clean
@cd ..; tar cvzf src/$(ROOT)_src.tar.gz \
src/Make* src/Package.sh src/MAKE src/*.cpp src/*.h src/STUBS \
$(patsubst %,src/%,$(PACKAGEUC)) $(patsubst %,src/%,$(PACKUSERUC)) \
--exclude=*/.svn
@cd STUBS; $(MAKE)
@echo "Created $(ROOT)_src.tar.gz"
# Make MPI STUBS library
stubs:
@cd STUBS; $(MAKE) clean; $(MAKE)
# Create Makefile.lib, Makefile.shlib, and Makefile.list
makelib:
@ -135,11 +121,27 @@ makelist:
@$(SHELL) Make.sh style
@$(SHELL) Make.sh Makefile.list
# Make MPI STUBS library
stubs:
@cd STUBS; $(MAKE) clean; $(MAKE)
# install LAMMPS shared lib and Python wrapper for Python usage
install-python:
@python ../python/install.py
# Create a tarball of src dir and packages
tar:
@cd STUBS; $(MAKE) clean
@cd ..; tar cvzf src/$(ROOT)_src.tar.gz \
src/Make* src/Package.sh src/MAKE src/*.cpp src/*.h src/STUBS \
$(patsubst %,src/%,$(PACKAGEUC)) $(patsubst %,src/%,$(PACKUSERUC)) \
--exclude=*/.svn
@cd STUBS; $(MAKE)
@echo "Created $(ROOT)_src.tar.gz"
# Package management
package:
@ -218,6 +220,7 @@ no-%:
# update = replace src files with newer package files
# overwrite = overwrite package files with newer src files
# diff = show differences between src and package files
# purge = delete obsolete and auto-generated package files
package-status ps:
@for p in $(PACKAGEUC); do $(SHELL) Package.sh $$p status; done
@ -238,3 +241,9 @@ package-diff pd:
@for p in $(PACKAGEUC); do $(SHELL) Package.sh $$p diff; done
@echo ''
@for p in $(PACKUSERUC); do $(SHELL) Package.sh $$p diff; done
package-purge: Purge.list
@echo 'Purging obsolete and auto-generated source files'
@for f in `grep -v '#' Purge.list` ; \
do test -f $$f && rm $$f && echo $$f || : ; \
done

276
src/SNAP/compute_sna_atom.cpp Normal file
View File

@ -0,0 +1,276 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "sna.h"
#include "string.h"
#include "stdlib.h"
#include "compute_sna_atom.h"
#include "atom.h"
#include "update.h"
#include "modify.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 "openmp_snap.h"
using namespace LAMMPS_NS;
ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
double rmin0, rfac0;
int twojmax, switchflag;
radelem = NULL;
wjelem = NULL;
int ntypes = atom->ntypes;
int nargmin = 6+2*ntypes;
if (narg < nargmin) error->all(FLERR,"Illegal compute sna/atom command");
// default values
diagonalstyle = 0;
rmin0 = 0.0;
switchflag = 1;
// offset by 1 to match up with types
memory->create(radelem,ntypes+1,"sna/atom:radelem");
memory->create(wjelem,ntypes+1,"sna/atom:wjelem");
rcutfac = atof(arg[3]);
rfac0 = atof(arg[4]);
twojmax = atoi(arg[5]);
for(int i = 0; i < ntypes; i++)
radelem[i+1] = atof(arg[6+i]);
for(int i = 0; i < ntypes; i++)
wjelem[i+1] = atof(arg[6+ntypes+i]);
// construct cutsq
double cut;
cutmax = 0.0;
memory->create(cutsq,ntypes+1,ntypes+1,"sna/atom:cutsq");
for(int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
if (cut > cutmax) cutmax = cut;
cutsq[i][i] = cut*cut;
for(int j = i+1; j <= ntypes; j++) {
cut = (radelem[i]+radelem[j])*rcutfac;
cutsq[i][j] = cutsq[j][i] = cut*cut;
}
}
// process optional args
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"diagonal") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
diagonalstyle = atoi(arg[iarg+1]);
if (diagonalstyle < 0 || diagonalstyle > 3)
error->all(FLERR,"Illegal compute sna/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
rmin0 = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"switchflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
switchflag = atoi(arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal compute sna/atom command");
}
// always unset use_shared_arrays
int use_shared_arrays = 0;
int nthreads = omp_get_max_threads();
snaptr = new SNA*[nthreads];
#pragma omp parallel shared(rfac0,twojmax,rmin0)
{
int tid = omp_get_thread_num();
snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,use_shared_arrays,
rmin0,switchflag);
}
ncoeff = snaptr[0]->ncoeff;
peratom_flag = 1;
size_peratom_cols = ncoeff;
nmax = 0;
njmax = 0;
sna = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSNAAtom::~ComputeSNAAtom()
{
memory->destroy(sna);
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(cutsq);
delete [] snaptr;
}
/* ---------------------------------------------------------------------- */
void ComputeSNAAtom::init()
{
if (force->pair == NULL)
error->all(FLERR,"Compute sna/atom requires a pair style be defined");
if (cutmax > force->pair->cutforce)
error->all(FLERR,"Compute sna/atom cutoff is longer than pairwise cutoff");
// 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,"sna/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute sna/atom");
#pragma omp parallel
{
int tid = omp_get_thread_num();
snaptr[tid]->init();
}
}
/* ---------------------------------------------------------------------- */
void ComputeSNAAtom::init_list(int id, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeSNAAtom::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow sna array if necessary
if (atom->nlocal > nmax) {
memory->destroy(sna);
nmax = atom->nmax;
memory->create(sna,nmax,size_peratom_cols,"sna/atom:sna");
array_atom = sna;
}
// invoke full neighbor list (will copy or build if necessary)
neighbor->build_one(list->index);
const int inum = list->inum;
const int* const ilist = list->ilist;
const int* const numneigh = list->numneigh;
int** const firstneigh = list->firstneigh;
int *type = atom->type;
// compute sna for each atom in group
// use full neighbor list to count atoms less than cutoff
double** const x = atom->x;
const int* const mask = atom->mask;
#pragma omp parallel for
for (int ii = 0; ii < inum; ii++) {
const int tid = omp_get_thread_num();
const int i = ilist[ii];
if (mask[i] & groupbit) {
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
const int itype = type[i];
const double radi = radelem[itype];
const int* const jlist = firstneigh[i];
const int jnum = numneigh[i];
// insure rij, inside, and typej are of size jnum
snaptr[tid]->grow_rij(jnum);
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
// typej = types of neighbors of I within cutoff
int ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j];
if (rsq < cutsq[itype][jtype] && rsq>1e-20) {
snaptr[tid]->rij[ninside][0] = delx;
snaptr[tid]->rij[ninside][1] = dely;
snaptr[tid]->rij[ninside][2] = delz;
snaptr[tid]->inside[ninside] = j;
snaptr[tid]->wj[ninside] = wjelem[jtype];
snaptr[tid]->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
ninside++;
}
}
snaptr[tid]->compute_ui(ninside);
snaptr[tid]->compute_zi();
snaptr[tid]->compute_bi();
snaptr[tid]->copy_bi2bvec();
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
sna[i][icoeff] = snaptr[tid]->bvec[icoeff];
} else {
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
sna[i][icoeff] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double ComputeSNAAtom::memory_usage()
{
double bytes = nmax*size_peratom_cols * sizeof(double);
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
bytes += snaptr[0]->memory_usage()*omp_get_max_threads();
return bytes;
}

53
src/SNAP/compute_sna_atom.h Normal file
View File

@ -0,0 +1,53 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(sna/atom,ComputeSNAAtom)
#else
#ifndef LMP_COMPUTE_SNA_ATOM_H
#define LMP_COMPUTE_SNA_ATOM_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSNAAtom : public Compute {
public:
ComputeSNAAtom(class LAMMPS *, int, char **);
~ComputeSNAAtom();
void init();
void init_list(int, class NeighList *);
void compute_peratom();
double memory_usage();
private:
int nmax, njmax, diagonalstyle;
int ncoeff;
double **cutsq;
class NeighList *list;
double **sna;
double rcutfac;
double *radelem;
double *wjelem;
class SNA** snaptr;
double cutmax;
};
}
#endif
#endif

329
src/SNAP/compute_snad_atom.cpp Normal file
View File

@ -0,0 +1,329 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "sna.h"
#include "string.h"
#include "stdlib.h"
#include "compute_snad_atom.h"
#include "atom.h"
#include "update.h"
#include "modify.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 "openmp_snap.h"
using namespace LAMMPS_NS;
ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
double rfac0, rmin0;
int twojmax, switchflag;
radelem = NULL;
wjelem = NULL;
int ntypes = atom->ntypes;
int nargmin = 6+2*ntypes;
if (narg < nargmin) error->all(FLERR,"Illegal compute snad/atom command");
// default values
diagonalstyle = 0;
rmin0 = 0;
switchflag = 1;
// process required arguments
memory->create(radelem,ntypes+1,"sna/atom:radelem"); // offset by 1 to match up with types
memory->create(wjelem,ntypes+1,"sna/atom:wjelem");
rcutfac = atof(arg[3]);
rfac0 = atof(arg[4]);
twojmax = atoi(arg[5]);
for(int i = 0; i < ntypes; i++)
radelem[i+1] = atof(arg[6+i]);
for(int i = 0; i < ntypes; i++)
wjelem[i+1] = atof(arg[6+ntypes+i]);
// construct cutsq
double cut;
memory->create(cutsq,ntypes+1,ntypes+1,"sna/atom:cutsq");
for(int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
cutsq[i][i] = cut*cut;
for(int j = i+1; j <= ntypes; j++) {
cut = (radelem[i]+radelem[j])*rcutfac;
cutsq[i][j] = cutsq[j][i] = cut*cut;
}
}
// process optional args
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"diagonal") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
diagonalstyle = atof(arg[iarg+1]);
if (diagonalstyle < 0 || diagonalstyle > 3)
error->all(FLERR,"Illegal compute snad/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
rmin0 = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"switchflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
switchflag = atoi(arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal compute snad/atom command");
}
// use_shared_arrays doesn't work with computes and should be disabled.
int use_shared_arrays = 0;
int nthreads = omp_get_max_threads();
snaptr = new SNA*[nthreads];
#pragma omp parallel shared(rfac0,twojmax,rmin0)
{
int tid = omp_get_thread_num();
snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,use_shared_arrays,
rmin0,switchflag);
}
ncoeff = snaptr[0]->ncoeff;
peratom_flag = 1;
size_peratom_cols = 3*ncoeff*atom->ntypes;
comm_reverse = size_peratom_cols;
nmax = 0;
njmax = 0;
snad = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSNADAtom::~ComputeSNADAtom()
{
memory->destroy(snad);
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(cutsq);
delete [] snaptr;
}
/* ---------------------------------------------------------------------- */
void ComputeSNADAtom::init()
{
if (force->pair == NULL)
error->all(FLERR,"Compute snad/atom requires a pair style be defined");
// TODO: Not sure what to do with this error check since cutoff radius is not
// a single number
//if (sqrt(cutsq) > force->pair->cutforce)
//error->all(FLERR,"Compute snad/atom cutoff is longer than pairwise cutoff");
// 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,"snad/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute snad/atom");
#pragma omp parallel
{
int tid = omp_get_thread_num();
snaptr[tid]->init();
}
}
/* ---------------------------------------------------------------------- */
void ComputeSNADAtom::init_list(int id, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeSNADAtom::compute_peratom()
{
int ntotal = atom->nlocal + atom->nghost;
invoked_peratom = update->ntimestep;
// grow snad array if necessary
if (ntotal > nmax) {
memory->destroy(snad);
nmax = atom->nmax;
memory->create(snad,nmax,size_peratom_cols,
"snad/atom:snad");
array_atom = snad;
}
// clear local array
for (int i = 0; i < ntotal; i++)
for (int icoeff = 0; icoeff < size_peratom_cols; icoeff++) {
snad[i][icoeff] = 0.0;
}
// invoke full neighbor list (will copy or build if necessary)
neighbor->build_one(list->index);
const int inum = list->inum;
const int* const ilist = list->ilist;
const int* const numneigh = list->numneigh;
int** const firstneigh = list->firstneigh;
int *type = atom->type;
// compute sna derivatives for each atom in group
// use full neighbor list to count atoms less than cutoff
double** const x = atom->x;
const int* const mask = atom->mask;
#pragma omp parallel for
for (int ii = 0; ii < inum; ii++) {
const int tid = omp_get_thread_num();
const int i = ilist[ii];
if (mask[i] & groupbit) {
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
const int itype = type[i];
const double radi = radelem[itype];
const int* const jlist = firstneigh[i];
const int jnum = numneigh[i];
const int typeoffset = 3*ncoeff*(atom->type[i]-1);
// insure rij, inside, and typej are of size jnum
snaptr[tid]->grow_rij(jnum);
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
// typej = types of neighbors of I within cutoff
// note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
int ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double delx = x[j][0] - xtmp;
const double dely = x[j][1] - ytmp;
const double delz = x[j][2] - ztmp;
const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
snaptr[tid]->rij[ninside][0] = delx;
snaptr[tid]->rij[ninside][1] = dely;
snaptr[tid]->rij[ninside][2] = delz;
snaptr[tid]->inside[ninside] = j;
snaptr[tid]->wj[ninside] = wjelem[jtype];
snaptr[tid]->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
ninside++;
}
}
snaptr[tid]->compute_ui(ninside);
snaptr[tid]->compute_zi();
for (int jj = 0; jj < ninside; jj++) {
const int j = snaptr[tid]->inside[jj];
snaptr[tid]->compute_duidrj(snaptr[tid]->rij[jj],
snaptr[tid]->wj[jj],
snaptr[tid]->rcutij[jj]);
snaptr[tid]->compute_dbidrj();
snaptr[tid]->copy_dbi2dbvec();
// Accumulate -dBi/dRi, -dBi/dRj
double *snadi = snad[i]+typeoffset;
double *snadj = snad[j]+typeoffset;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
snadi[icoeff] += snaptr[tid]->dbvec[icoeff][0];
snadi[icoeff+ncoeff] += snaptr[tid]->dbvec[icoeff][1];
snadi[icoeff+2*ncoeff] += snaptr[tid]->dbvec[icoeff][2];
snadj[icoeff] -= snaptr[tid]->dbvec[icoeff][0];
snadj[icoeff+ncoeff] -= snaptr[tid]->dbvec[icoeff][1];
snadj[icoeff+2*ncoeff] -= snaptr[tid]->dbvec[icoeff][2];
}
}
}
}
// communicate snad contributions between neighbor procs
comm->reverse_comm_compute(this);
}
/* ---------------------------------------------------------------------- */
int ComputeSNADAtom::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last,icoeff;
m = 0;
last = first + n;
for (i = first; i < last; i++)
for (icoeff = 0; icoeff < size_peratom_cols; icoeff++)
buf[m++] = snad[i][icoeff];
return comm_reverse;
}
/* ---------------------------------------------------------------------- */
void ComputeSNADAtom::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m,icoeff;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
for (icoeff = 0; icoeff < size_peratom_cols; icoeff++)
snad[j][icoeff] += buf[m++];
}
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double ComputeSNADAtom::memory_usage()
{
double bytes = nmax*size_peratom_cols * sizeof(double);
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
bytes += ncoeff*3;
bytes += snaptr[0]->memory_usage()*omp_get_max_threads();
return bytes;
}

54
src/SNAP/compute_snad_atom.h Normal file
View File

@ -0,0 +1,54 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(snad/atom,ComputeSNADAtom)
#else
#ifndef LMP_COMPUTE_SNAD_ATOM_H
#define LMP_COMPUTE_SNAD_ATOM_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSNADAtom : public Compute {
public:
ComputeSNADAtom(class LAMMPS *, int, char **);
~ComputeSNADAtom();
void init();
void init_list(int, class NeighList *);
void compute_peratom();
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
double memory_usage();
private:
int nmax, njmax, diagonalstyle;
int ncoeff;
double **cutsq;
class NeighList *list;
double **snad;
double rcutfac;
double *radelem;
double *wjelem;
class SNA** snaptr;
};
}
#endif
#endif

339
src/SNAP/compute_snav_atom.cpp Normal file
View File

@ -0,0 +1,339 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "sna.h"
#include "string.h"
#include "stdlib.h"
#include "compute_snav_atom.h"
#include "atom.h"
#include "update.h"
#include "modify.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 "openmp_snap.h"
using namespace LAMMPS_NS;
ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
double rfac0, rmin0;
int twojmax, switchflag;
radelem = NULL;
wjelem = NULL;
nvirial = 6;
int ntypes = atom->ntypes;
int nargmin = 6+2*ntypes;
if (narg < nargmin) error->all(FLERR,"Illegal compute snav/atom command");
// default values
diagonalstyle = 0;
rmin0 = 0;
switchflag = 1;
// process required arguments
memory->create(radelem,ntypes+1,"sna/atom:radelem"); // offset by 1 to match up with types
memory->create(wjelem,ntypes+1,"sna/atom:wjelem");
rcutfac = atof(arg[3]);
rfac0 = atof(arg[4]);
twojmax = atoi(arg[5]);
for(int i = 0; i < ntypes; i++)
radelem[i+1] = atof(arg[6+i]);
for(int i = 0; i < ntypes; i++)
wjelem[i+1] = atof(arg[6+ntypes+i]);
// construct cutsq
double cut;
memory->create(cutsq,ntypes+1,ntypes+1,"sna/atom:cutsq");
for(int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
cutsq[i][i] = cut*cut;
for(int j = i+1; j <= ntypes; j++) {
cut = (radelem[i]+radelem[j])*rcutfac;
cutsq[i][j] = cutsq[j][i] = cut*cut;
}
}
// process optional args
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"diagonal") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
diagonalstyle = atof(arg[iarg+1]);
if (diagonalstyle < 0 || diagonalstyle > 3)
error->all(FLERR,"Illegal compute snav/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
rmin0 = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"switchflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
switchflag = atoi(arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal compute snav/atom command");
}
// use_shared_arrays doesn't work with computes and should be disabled.
int use_shared_arrays = 0;
int nthreads = omp_get_max_threads();
snaptr = new SNA*[nthreads];
#pragma omp parallel shared(rfac0,twojmax,rmin0,switchflag)
{
int tid = omp_get_thread_num();
snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,use_shared_arrays,
rmin0,switchflag);
}
ncoeff = snaptr[0]->ncoeff;
peratom_flag = 1;
size_peratom_cols = nvirial*ncoeff*atom->ntypes;
comm_reverse = size_peratom_cols;
nmax = 0;
njmax = 0;
snav = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSNAVAtom::~ComputeSNAVAtom()
{
memory->destroy(snav);
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(cutsq);
delete [] snaptr;
}
/* ---------------------------------------------------------------------- */
void ComputeSNAVAtom::init()
{
if (force->pair == NULL)
error->all(FLERR,"Compute snav/atom requires a pair style be defined");
// TODO: Not sure what to do with this error check since cutoff radius is not
// a single number
//if (sqrt(cutsq) > force->pair->cutforce)
// error->all(FLERR,"Compute snav/atom cutoff is longer than pairwise cutoff");
// 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,"snav/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute snav/atom");
#pragma omp parallel
{
int tid = omp_get_thread_num();
snaptr[tid]->init();
}
}
/* ---------------------------------------------------------------------- */
void ComputeSNAVAtom::init_list(int id, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeSNAVAtom::compute_peratom()
{
int ntotal = atom->nlocal + atom->nghost;
invoked_peratom = update->ntimestep;
// grow snav array if necessary
if (ntotal > nmax) {
memory->destroy(snav);
nmax = atom->nmax;
memory->create(snav,nmax,size_peratom_cols,
"snav/atom:snav");
array_atom = snav;
}
// clear local array
for (int i = 0; i < ntotal; i++)
for (int icoeff = 0; icoeff < size_peratom_cols; icoeff++) {
snav[i][icoeff] = 0.0;
}
// invoke full neighbor list (will copy or build if necessary)
neighbor->build_one(list->index);
const int inum = list->inum;
const int* const ilist = list->ilist;
const int* const numneigh = list->numneigh;
int** const firstneigh = list->firstneigh;
int *type = atom->type;
// compute sna derivatives for each atom in group
// use full neighbor list to count atoms less than cutoff
double** const x = atom->x;
const int* const mask = atom->mask;
#pragma omp parallel for
for (int ii = 0; ii < inum; ii++) {
const int tid = omp_get_thread_num();
const int i = ilist[ii];
if (mask[i] & groupbit) {
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
const int itype = type[i];
const double radi = radelem[itype];
const int* const jlist = firstneigh[i];
const int jnum = numneigh[i];
const int typeoffset = nvirial*ncoeff*(atom->type[i]-1);
// insure rij, inside, and typej are of size jnum
snaptr[tid]->grow_rij(jnum);
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
// typej = types of neighbors of I within cutoff
// note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
int ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double delx = x[j][0] - xtmp;
const double dely = x[j][1] - ytmp;
const double delz = x[j][2] - ztmp;
const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
snaptr[tid]->rij[ninside][0] = delx;
snaptr[tid]->rij[ninside][1] = dely;
snaptr[tid]->rij[ninside][2] = delz;
snaptr[tid]->inside[ninside] = j;
snaptr[tid]->wj[ninside] = wjelem[jtype];
snaptr[tid]->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
ninside++;
}
}
snaptr[tid]->compute_ui(ninside);
snaptr[tid]->compute_zi();
for (int jj = 0; jj < ninside; jj++) {
const int j = snaptr[tid]->inside[jj];
snaptr[tid]->compute_duidrj(snaptr[tid]->rij[jj],
snaptr[tid]->wj[jj],
snaptr[tid]->rcutij[jj]);
snaptr[tid]->compute_dbidrj();
snaptr[tid]->copy_dbi2dbvec();
// Accumulate -dBi/dRi*Ri, -dBi/dRj*Rj
double *snavi = snav[i]+typeoffset;
double *snavj = snav[j]+typeoffset;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
snavi[icoeff] += snaptr[tid]->dbvec[icoeff][0]*xtmp;
snavi[icoeff+ncoeff] += snaptr[tid]->dbvec[icoeff][1]*ytmp;
snavi[icoeff+2*ncoeff] += snaptr[tid]->dbvec[icoeff][2]*ztmp;
snavi[icoeff+3*ncoeff] += snaptr[tid]->dbvec[icoeff][1]*ztmp;
snavi[icoeff+4*ncoeff] += snaptr[tid]->dbvec[icoeff][0]*ztmp;
snavi[icoeff+5*ncoeff] += snaptr[tid]->dbvec[icoeff][0]*ytmp;
snavj[icoeff] -= snaptr[tid]->dbvec[icoeff][0]*x[j][0];
snavj[icoeff+ncoeff] -= snaptr[tid]->dbvec[icoeff][1]*x[j][1];
snavj[icoeff+2*ncoeff] -= snaptr[tid]->dbvec[icoeff][2]*x[j][2];
snavj[icoeff+3*ncoeff] -= snaptr[tid]->dbvec[icoeff][1]*x[j][2];
snavj[icoeff+4*ncoeff] -= snaptr[tid]->dbvec[icoeff][0]*x[j][2];
snavj[icoeff+5*ncoeff] -= snaptr[tid]->dbvec[icoeff][0]*x[j][1];
}
}
}
}
// communicate snav contributions between neighbor procs
comm->reverse_comm_compute(this);
}
/* ---------------------------------------------------------------------- */
int ComputeSNAVAtom::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last,icoeff;
m = 0;
last = first + n;
for (i = first; i < last; i++)
for (icoeff = 0; icoeff < size_peratom_cols; icoeff++)
buf[m++] = snav[i][icoeff];
return comm_reverse;
}
/* ---------------------------------------------------------------------- */
void ComputeSNAVAtom::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m,icoeff;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
for (icoeff = 0; icoeff < size_peratom_cols; icoeff++)
snav[j][icoeff] += buf[m++];
}
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double ComputeSNAVAtom::memory_usage()
{
double bytes = nmax*size_peratom_cols * sizeof(double);
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
bytes += ncoeff*nvirial;
bytes += snaptr[0]->memory_usage()*omp_get_max_threads();
return bytes;
}

55
src/SNAP/compute_snav_atom.h Normal file
View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(snav/atom,ComputeSNAVAtom)
#else
#ifndef LMP_COMPUTE_SNAV_ATOM_H
#define LMP_COMPUTE_SNAV_ATOM_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSNAVAtom : public Compute {
public:
ComputeSNAVAtom(class LAMMPS *, int, char **);
~ComputeSNAVAtom();
void init();
void init_list(int, class NeighList *);
void compute_peratom();
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
double memory_usage();
private:
int nmax, njmax, diagonalstyle;
int ncoeff,nvirial;
double **cutsq;
class NeighList *list;
double **snav;
double rcutfac;
double *radelem;
double *wjelem;
class SNA** snaptr;
};
}
#endif
#endif

11
src/SNAP/openmp_snap.h Normal file
View File

@ -0,0 +1,11 @@
#ifdef _OPENMP
#include <omp.h>
#else
enum omp_sched_t {omp_sched_static, omp_sched_dynamic, omp_sched_guided, omp_sched_auto};
inline int omp_get_thread_num() { return 0;}
inline int omp_get_max_threads() { return 1;}
inline int omp_set_num_threads(int num_threads) {return 1;}
/* inline int __sync_fetch_and_add(int* ptr, int value) {int tmp = *ptr; ptr[0]+=value; return tmp;} */
inline void omp_set_schedule(omp_sched_t schedule,int modifier=1) {}
inline int omp_in_parallel() {return 0;}
#endif

1725
src/SNAP/pair_snap.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

107
src/SNAP/pair_snap.h Normal file
View File

@ -0,0 +1,107 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(snap,PairSNAP)
#else
#ifndef LMP_PAIR_SNAP_H
#define LMP_PAIR_SNAP_H
#include "pair.h"
namespace LAMMPS_NS {
class PairSNAP : public Pair {
public:
PairSNAP(class LAMMPS *);
~PairSNAP();
void compute(int, int);
void compute_regular(int, int);
void compute_optimized(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
double memory_usage();
protected:
int ncoeff;
double **bvec, ***dbvec;
class SNA** sna;
int nmax;
int nthreads;
void allocate();
void read_files(char *, char *);
inline int equal(double* x,double* y);
inline double dist2(double* x,double* y);
double extra_cutoff();
void load_balance();
void set_sna_to_shared(int snaid,int i);
void build_per_atom_arrays();
int schedule_user;
double schedule_time_guided;
double schedule_time_dynamic;
int ncalls_neigh;
int do_load_balance;
int ilistmask_max;
int* ilistmask;
int ghostinum;
int ghostilist_max;
int* ghostilist;
int ghostnumneigh_max;
int* ghostnumneigh;
int* ghostneighs;
int* ghostfirstneigh;
int ghostneighs_total;
int ghostneighs_max;
int use_optimized;
int use_shared_arrays;
int i_max;
int i_neighmax;
int i_numpairs;
int **i_pairs;
double ***i_rij;
int **i_inside;
double **i_wj;
double **i_rcutij;
int *i_ninside;
double ****i_uarraytot_r, ****i_uarraytot_i;
double ******i_zarray_r, ******i_zarray_i;
// timespec starttime, endtime;
double timers[4];
double gamma;
double rcutmax; // max cutoff for all elements
int nelements; // # of unique elements
char **elements; // names of unique elements
double *radelem; // element radii
double *wjelem; // elements weights
double **coeffelem; // element bispectrum coefficients
int *map; // mapping from atom types to elements
int twojmax, diagonalstyle, switchflag;
double rcutfac, rfac0, rmin0, wj1, wj2;
int rcutfacflag, twojmaxflag; // flags for required parameters
int gammaoneflag; // 1 if parameter gamma is 1
};
}
#endif
#endif

2082
src/SNAP/sna.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

135
src/SNAP/sna.h Normal file
View File

@ -0,0 +1,135 @@
/* ----------------------------------------------------------------------
Authors: Aidan Thompson and Christian Trott, Sandia National Labs, 2012
Property of Sandia National Labs: Not for External Distribution
------------------------------------------------------------------------- */
#ifndef LMP_SNA_H
#define LMP_SNA_H
#include <complex>
#include "pointers.h"
#include <ctime>
namespace LAMMPS_NS {
struct SNA_LOOPINDICES {
int j1, j2, j, ma, mb, ma1, ma2, mb1, mb2;
};
struct SNA_LOOPINDICES_J {
int j1, j2, j;
};
class SNA : protected Pointers {
public:
SNA(LAMMPS*, double, int, int, int, double, int);
SNA(LAMMPS* lmp) : Pointers(lmp) {};
~SNA();
void build_indexlist();
void init();
void test();
double memory_usage();
int ncoeff;
// functions for bispectrum coefficients
void compute_ui(int);
void compute_ui_omp(int, int);
void compute_zi();
void compute_zi_omp(int);
void compute_bi();
void copy_bi2bvec();
// functions for derivatives
void compute_duidrj(double*, double, double);
void compute_dbidrj();
void compute_dbidrj_nonsymm();
void copy_dbi2dbvec();
double compute_sfac(double, double);
double compute_dsfac(double, double);
double* timers;
timespec starttime, endtime;
int print;
int counter;
//per sna class instance for OMP use
double* bvec, ** dbvec;
double** rij;
int* inside;
double* wj;
double* rcutij;
int nmax;
void grow_rij(int);
int twojmax, diagonalstyle;
double*** uarraytot_r, *** uarraytot_i;
double***** zarray_r, ***** zarray_i;
double*** uarraytot_r_b, *** uarraytot_i_b;
double***** zarray_r_b, ***** zarray_i_b;
double*** uarray_r, *** uarray_i;
private:
double rmin0, rfac0;
//use indexlist instead of loops, constructor generates these
SNA_LOOPINDICES* idx;
int idx_max;
SNA_LOOPINDICES_J* idxj;
int idxj_max;
// data for bispectrum coefficients
double***** cgarray;
double** rootpqarray;
double*** barray;
// derivatives of data
double**** duarray_r, **** duarray_i;
double**** dbarray;
void create_twojmax_arrays();
void destroy_twojmax_arrays();
void init_clebsch_gordan();
void init_rootpqarray();
void jtostr(char*, int);
void mtostr(char*, int, int);
void print_clebsch_gordan(FILE*);
void zero_uarraytot();
void addself_uarraytot(double);
void add_uarraytot(double, double, double);
void add_uarraytot_omp(double, double, double);
void compute_uarray(double, double, double,
double, double);
void compute_uarray_omp(double, double, double,
double, double, int);
double factorial(int);
double deltacg(int, int, int);
int compute_ncoeff();
void compute_duarray(double, double, double,
double, double, double, double, double);
// if number of atoms are small use per atom arrays
// for twojmax arrays, rij, inside, bvec
// this will increase the memory footprint considerably,
// but allows parallel filling and reuse of these arrays
int use_shared_arrays;
// Sets the style for the switching function
// 0 = none
// 1 = cosine
int switch_flag;
// Self-weight
double wself;
};
}
#endif