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

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

Browse Source
This commit is contained in:
sjplimp
2016-04-22 15:55:54 +00:00
parent a4af037d1d
commit f70a764d5c
13 changed files with 495 additions and 45 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
src/GPU/pair_coul_long_gpu.cpp
View File

@ -144,7 +144,14 @@ void PairCoulLongGPU::init_style()
if (force->newton_pair)
error->all(FLERR,"Cannot use newton pair with coul/long/gpu pair style");
// Repeat cutsq calculation because done after call to init_style
// Call init_one calculation make sure scale is correct
for (int i = 1; i <= atom->ntypes; i++) {
for (int j = i; j <= atom->ntypes; j++) {
if (setflag[i][j] != 0 || (setflag[i][i] != 0 && setflag[j][j] != 0)) {
double cut = init_one(i,j);
}
}
}
double cell_size = cut_coul + neighbor->skin;
cut_coulsq = cut_coul * cut_coul;

64
src/MOLECULE/bond_table.cpp
View File

@ -29,9 +29,10 @@
using namespace LAMMPS_NS;
enum{LINEAR,SPLINE};
enum{NONE,LINEAR,SPLINE};
#define MAXLINE 1024
#define BIGNUM 1.0e300
/* ---------------------------------------------------------------------- */
@ -134,6 +135,7 @@ void BondTable::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,"Illegal bond_style command");
tabstyle = NONE;
if (strcmp(arg[0],"linear") == 0) tabstyle = LINEAR;
else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
else error->all(FLERR,"Unknown table style in bond style table");
@ -289,6 +291,7 @@ void BondTable::free_table(Table *tb)
void BondTable::read_table(Table *tb, char *file, char *keyword)
{
char line[MAXLINE];
double emin = BIGNUM;
// open file
@ -326,14 +329,65 @@ void BondTable::read_table(Table *tb, char *file, char *keyword)
// read r,e,f table values from file
int itmp;
int cerror = 0;
int r0idx = -1;
fgets(line,MAXLINE,fp);
for (int i = 0; i < tb->ninput; i++) {
fgets(line,MAXLINE,fp);
sscanf(line,"%d %lg %lg %lg",
&itmp,&tb->rfile[i],&tb->efile[i],&tb->ffile[i]);
if (NULL == fgets(line,MAXLINE,fp))
error->one(FLERR,"Premature end of file in bond table");
if (4 != sscanf(line,"%d %lg %lg %lg",
&itmp,&tb->rfile[i],&tb->efile[i],&tb->ffile[i])) ++cerror;
if (tb->efile[i] < emin) {
emin = tb->efile[i];
r0idx = i;
}
}
fclose(fp);
// infer r0 from minimum of potential, if not given explicitly
if ((tb->r0 == 0.0) && (r0idx >= 0)) tb->r0 = tb->rfile[r0idx];
// warn if force != dE/dr at any point that is not an inflection point
// check via secant approximation to dE/dr
// skip two end points since do not have surrounding secants
// inflection point is where curvature changes sign
double r,e,f,rprev,rnext,eprev,enext,fleft,fright;
int ferror = 0;
for (int i = 1; i < tb->ninput-1; i++) {
r = tb->rfile[i];
rprev = tb->rfile[i-1];
rnext = tb->rfile[i+1];
e = tb->efile[i];
eprev = tb->efile[i-1];
enext = tb->efile[i+1];
f = tb->ffile[i];
fleft = - (e-eprev) / (r-rprev);
fright = - (enext-e) / (rnext-r);
if (f < fleft && f < fright) ferror++;
if (f > fleft && f > fright) ferror++;
//printf("Values %d: %g %g %g\n",i,r,e,f);
//printf(" secant %d %d %g: %g %g %g\n",i,ferror,r,fleft,fright,f);
}
fclose(fp);
if (ferror) {
char str[128];
sprintf(str,"%d of %d force values in table are inconsistent with -dE/dr.\n"
" Should only be flagged at inflection points",ferror,tb->ninput);
error->warning(FLERR,str);
}
// warn if data was read incompletely, e.g. columns were missing
if (cerror) {
char str[128];
sprintf(str,"%d of %d lines in table were incomplete or could not be"
" parsed completely",cerror,tb->ninput);
error->warning(FLERR,str);
}
}
/* ----------------------------------------------------------------------

5
src/USER-FEP/README
View File

@ -11,3 +11,8 @@ There are auxiliary tools for using this package in tools/fep.
The person who created this package is Agilio Padua at Université
Blaise Pascal Clermont-Ferrand (agilio.padua at univ-bpclermont.fr)
Contact him directly if you have questions.
Pair style morse/soft was contributed by Stefan Paquay, s.paquay at tue.nl
Applied Physics/Theory of Polymers and Soft Matter,
Eindhoven University of Technology (TU/e), The Netherlands
Contact him in case of problems with this pair style.

295
src/compute_dipole_chunk.cpp Normal file
View File

@ -0,0 +1,295 @@
/* ----------------------------------------------------------------------
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 "string.h"
#include "compute_dipole_chunk.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#include "math_special.h"
using namespace LAMMPS_NS;
using namespace MathSpecial;
enum { MASSCENTER, GEOMCENTER };
/* ---------------------------------------------------------------------- */
ComputeDipoleChunk::ComputeDipoleChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if ((narg != 4) && (narg != 5)) error->all(FLERR,"Illegal compute dipole/chunk command");
array_flag = 1;
size_array_cols = 4;
size_array_rows = 0;
size_array_rows_variable = 1;
extarray = 0;
// ID of compute chunk/atom
int n = strlen(arg[3]) + 1;
idchunk = new char[n];
strcpy(idchunk,arg[3]);
usecenter = MASSCENTER;
if (narg == 5) {
if (strncmp(arg[4],"geom",4) == 0) usecenter = GEOMCENTER;
else if (strcmp(arg[4],"mass") == 0) usecenter = MASSCENTER;
else error->all(FLERR,"Illegal compute dipole/chunk command");
}
init();
// chunk-based data
nchunk = 1;
maxchunk = 0;
massproc = masstotal = NULL;
chrgproc = chrgtotal = NULL;
com = comall = NULL;
dipole = dipoleall = NULL;
allocate();
}
/* ---------------------------------------------------------------------- */
ComputeDipoleChunk::~ComputeDipoleChunk()
{
delete [] idchunk;
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(chrgproc);
memory->destroy(chrgtotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(dipole);
memory->destroy(dipoleall);
}
/* ---------------------------------------------------------------------- */
void ComputeDipoleChunk::init()
{
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for "
"compute dipole/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute dipole/chunk does not use chunk/atom compute");
}
/* ---------------------------------------------------------------------- */
void ComputeDipoleChunk::compute_array()
{
int i,index;
double massone;
double unwrap[3];
invoked_array = update->ntimestep;
// compute chunk/atom assigns atoms to chunk IDs
// extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
nchunk = cchunk->setup_chunks();
cchunk->compute_ichunk();
int *ichunk = cchunk->ichunk;
if (nchunk > maxchunk) allocate();
size_array_rows = nchunk;
// zero local per-chunk values
for (int i = 0; i < nchunk; i++) {
massproc[i] = chrgproc[i] = 0.0;
com[i][0] = com[i][1] = com[i][2] = 0.0;
dipole[i][0] = dipole[i][1] = dipole[i][2] = dipole[i][3] = 0.0;
}
// compute COM for each chunk
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
double *q = atom->q;
double **mu = atom->mu;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
if (usecenter == MASSCENTER) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
} else massone = 1.0; // usecenter == GEOMCENTER
domain->unmap(x[i],image[i],unwrap);
massproc[index] += massone;
if (atom->q_flag) chrgproc[index] += atom->q[i];
com[index][0] += unwrap[0] * massone;
com[index][1] += unwrap[1] * massone;
com[index][2] += unwrap[2] * massone;
}
MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(chrgproc,chrgtotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nchunk; i++) {
if (masstotal[i] == 0.0) masstotal[i] = 1.0;
comall[i][0] /= masstotal[i];
comall[i][1] /= masstotal[i];
comall[i][2] /= masstotal[i];
}
// compute dipole for each chunk
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
domain->unmap(x[i],image[i],unwrap);
if (atom->q_flag) {
dipole[index][0] += q[i]*unwrap[0];
dipole[index][1] += q[i]*unwrap[1];
dipole[index][2] += q[i]*unwrap[2];
}
if (atom->mu_flag) {
dipole[index][0] += mu[i][0];
dipole[index][1] += mu[i][1];
dipole[index][2] += mu[i][2];
}
}
}
MPI_Allreduce(&dipole[0][0],&dipoleall[0][0],4*nchunk,
MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < nchunk; i++) {
// correct for position dependence with charged chunks
dipoleall[i][0] -= chrgtotal[i]*comall[i][0];
dipoleall[i][1] -= chrgtotal[i]*comall[i][1];
dipoleall[i][2] -= chrgtotal[i]*comall[i][2];
// compute total dipole moment
dipoleall[i][3] = sqrt(square(dipoleall[i][0])
+ square(dipoleall[i][1])
+ square(dipoleall[i][2]));
}
}
/* ----------------------------------------------------------------------
lock methods: called by fix ave/time
these methods insure vector/array size is locked for Nfreq epoch
by passing lock info along to compute chunk/atom
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
increment lock counter
------------------------------------------------------------------------- */
void ComputeDipoleChunk::lock_enable()
{
cchunk->lockcount++;
}
/* ----------------------------------------------------------------------
decrement lock counter in compute chunk/atom, it if still exists
------------------------------------------------------------------------- */
void ComputeDipoleChunk::lock_disable()
{
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
cchunk->lockcount--;
}
}
/* ----------------------------------------------------------------------
calculate and return # of chunks = length of vector/array
------------------------------------------------------------------------- */
int ComputeDipoleChunk::lock_length()
{
nchunk = cchunk->setup_chunks();
return nchunk;
}
/* ----------------------------------------------------------------------
set the lock from startstep to stopstep
------------------------------------------------------------------------- */
void ComputeDipoleChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
cchunk->lock(fixptr,startstep,stopstep);
}
/* ----------------------------------------------------------------------
unset the lock
------------------------------------------------------------------------- */
void ComputeDipoleChunk::unlock(Fix *fixptr)
{
cchunk->unlock(fixptr);
}
/* ----------------------------------------------------------------------
free and reallocate per-chunk arrays
------------------------------------------------------------------------- */
void ComputeDipoleChunk::allocate()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(chrgproc);
memory->destroy(chrgtotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(dipole);
memory->destroy(dipoleall);
maxchunk = nchunk;
memory->create(massproc,maxchunk,"dipole/chunk:massproc");
memory->create(masstotal,maxchunk,"dipole/chunk:masstotal");
memory->create(chrgproc,maxchunk,"dipole/chunk:chrgproc");
memory->create(chrgtotal,maxchunk,"dipole/chunk:chrgtotal");
memory->create(com,maxchunk,3,"dipole/chunk:com");
memory->create(comall,maxchunk,3,"dipole/chunk:comall");
memory->create(dipole,maxchunk,4,"dipole/chunk:dipole");
memory->create(dipoleall,maxchunk,4,"dipole/chunk:dipoleall");
array = dipoleall;
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeDipoleChunk::memory_usage()
{
double bytes = (bigint) maxchunk * 2 * sizeof(double);
bytes += (bigint) maxchunk * 2*3 * sizeof(double);
bytes += (bigint) maxchunk * 2*4 * sizeof(double);
return bytes;
}

77
src/compute_dipole_chunk.h Normal file
View File

@ -0,0 +1,77 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(dipole/chunk,ComputeDipoleChunk)
#else
#ifndef LMP_COMPUTE_DIPOLE_CHUNK_H
#define LMP_COMPUTE_DIPOLE_CHUNK_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeDipoleChunk : public Compute {
public:
ComputeDipoleChunk(class LAMMPS *, int, char **);
~ComputeDipoleChunk();
void init();
void compute_array();
void lock_enable();
void lock_disable();
int lock_length();
void lock(class Fix *, bigint, bigint);
void unlock(class Fix *);
double memory_usage();
private:
int nchunk,maxchunk;
char *idchunk;
class ComputeChunkAtom *cchunk;
double *massproc,*masstotal;
double *chrgproc,*chrgtotal;
double **com,**comall;
double **dipole,**dipoleall;
int usecenter;
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Chunk/atom compute does not exist for compute dipole/chunk
Self-explanatory.
E: Compute dipole/chunk does not use chunk/atom compute
The style of the specified compute is not chunk/atom.
*/

8
src/info.cpp
View File

@ -27,6 +27,7 @@
#include "fix.h"
#include "force.h"
#include "pair.h"
#include "pair_hybrid.h"
#include "group.h"
#include "input.h"
#include "modify.h"
@ -281,6 +282,13 @@ void Info::command(int narg, char **arg)
fprintf(out,"Atoms = " BIGINT_FORMAT ", types = %d, style = %s\n",
atom->natoms, atom->ntypes, force->pair_style);
if (force->pair && strstr(force->pair_style,"hybrid")) {
PairHybrid *hybrid = (PairHybrid *)force->pair;
fprintf(out,"Hybrid sub-styles:");
for (int i=0; i < hybrid->nstyles; ++i)
fprintf(out," %s", hybrid->keywords[i]);
fputc('\n',out);
}
if (atom->molecular > 0) {
const char *msg;
msg = force->bond_style ? force->bond_style : "none";

6
src/library.cpp
View File

@ -417,7 +417,8 @@ void lammps_gather_atoms(void *ptr, char *name,
int flag = 0;
if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;
if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
if (flag && lmp->comm->me == 0) {
if (flag) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,"Library error in lammps_gather_atoms");
return;
}
@ -506,7 +507,8 @@ void lammps_scatter_atoms(void *ptr, char *name,
if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;
if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
if (lmp->atom->map_style == 0) flag = 1;
if (flag && lmp->comm->me == 0) {
if (flag) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms");
return;
}

22
src/pair.cpp
View File

@ -38,11 +38,11 @@
#include "suffix.h"
#include "atom_masks.h"
#include "memory.h"
#include "math_const.h"
#include "error.h"
using namespace LAMMPS_NS;
#define EWALD_F 1.12837917
using namespace MathConst;
enum{NONE,RLINEAR,RSQ,BMP};
@ -56,8 +56,6 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp)
{
instance_me = instance_total++;
THIRD = 1.0/3.0;
eng_vdwl = eng_coul = 0.0;
comm_forward = comm_reverse = comm_reverse_off = 0;
@ -385,7 +383,7 @@ void Pair::init_tables(double cut_coul, double *cut_respa)
ftable[i] = qqrd2e/r * fgamma;
etable[i] = qqrd2e/r * egamma;
} else {
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
etable[i] = qqrd2e/r * derfc;
}
} else {
@ -397,9 +395,9 @@ void Pair::init_tables(double cut_coul, double *cut_respa)
etable[i] = qqrd2e/r * egamma;
vtable[i] = qqrd2e/r * fgamma;
} else {
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2 - 1.0);
etable[i] = qqrd2e/r * derfc;
vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
vtable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
}
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
@ -408,7 +406,7 @@ void Pair::init_tables(double cut_coul, double *cut_respa)
ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
} else {
if (msmflag) ftable[i] = qqrd2e/r * fgamma;
else ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
else ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
ctable[i] = qqrd2e/r;
}
}
@ -481,7 +479,7 @@ void Pair::init_tables(double cut_coul, double *cut_respa)
f_tmp = qqrd2e/r * fgamma;
e_tmp = qqrd2e/r * egamma;
} else {
f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
e_tmp = qqrd2e/r * derfc;
}
} else {
@ -492,9 +490,9 @@ void Pair::init_tables(double cut_coul, double *cut_respa)
e_tmp = qqrd2e/r * egamma;
v_tmp = qqrd2e/r * fgamma;
} else {
f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2 - 1.0);
e_tmp = qqrd2e/r * derfc;
v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
v_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
}
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
@ -503,7 +501,7 @@ void Pair::init_tables(double cut_coul, double *cut_respa)
c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
} else {
if (msmflag) f_tmp = qqrd2e/r * fgamma;
else f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
else f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
c_tmp = qqrd2e/r;
}
}

2
src/pair.h
View File

@ -219,8 +219,6 @@ class Pair : protected Pointers {
typedef union {int i; float f;} union_int_float_t;
double THIRD;
int vflag_fdotr;
int maxeatom,maxvatom;

1
src/pair_hybrid.h
View File

@ -31,6 +31,7 @@ class PairHybrid : public Pair {
friend class FixOMP;
friend class Force;
friend class Respa;
friend class Info;
public:
PairHybrid(class LAMMPS *);
virtual ~PairHybrid();

1
src/random_park.h
View File

@ -19,7 +19,6 @@
namespace LAMMPS_NS {
class RanPark : protected Pointers {
friend class Set;
public:
RanPark(class LAMMPS *, int);
double uniform();

5
src/write_restart.cpp
View File

@ -111,6 +111,7 @@ void WriteRestart::command(int narg, char **arg)
// init entire system since comm->exchange is done
// comm::init needs neighbor::init needs pair::init needs kspace::init, etc
if (noinit == 0) {
if (comm->me == 0 && screen)
fprintf(screen,"System init for write_restart ...\n");
lmp->init();
@ -133,6 +134,7 @@ void WriteRestart::command(int narg, char **arg)
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
}
// write single restart file
@ -220,6 +222,9 @@ void WriteRestart::multiproc_options(int multiproc_caller, int mpiioflag_caller,
else filewriter = 0;
iarg += 2;
} else if (strcmp(arg[iarg],"noinit") == 0) {
noinit = 1;
iarg++;
} else error->all(FLERR,"Illegal write_restart command");
}
}

1
src/write_restart.h
View File

@ -36,6 +36,7 @@ class WriteRestart : protected Pointers {
int me,nprocs;
FILE *fp;
bigint natoms; // natoms (sum of nlocal) to write into file
int noinit;
int multiproc; // 0 = proc 0 writes for all
// else # of procs writing files