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

Merge branch 'openmp-step8' into openmp-master

Browse Source
This commit is contained in:
Axel Kohlmeyer
2011-10-26 14:19:18 -04:00
parent 599c7027d4 3c110b35da
commit 793b6ddcaa
7 changed files with 430 additions and 537 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
doc/suffix.txt
View File

@ -37,9 +37,12 @@ USER-CUDA package.
These are the variants these packages provide:
OPT = a handful of pair styles, cache-optimized for faster CPU performance
USER-OMP = a collection of pair, dihedral and fix styles with support for OpenMP multi-threading
GPU = a handful of pair styles and the PPPM kspace_style, optimized to run on one or more GPUs or multicore CPU/GPU nodes
USER-CUDA = a collection of atom, pair, fix, compute, and intergrate styles, optimized to run on one or more NVIDIA GPUs :ul
USER-OMP = a collection of pair, bond, angle, dihedral, improper, kspace,
compute, and fix styles with support for OpenMP multi-threading
GPU = a handful of pair styles and the PPPM kspace_style, optimized to run
on one or more GPUs or multicore CPU/GPU nodes
USER-CUDA = a collection of atom, pair, fix, compute, and intergrate styles,
optimized to run on one or more NVIDIA GPUs :ul
As an example, all of the packages provide a "pair_style
lj/cut"_pair_lj.html variant, with style names lj/cut/opt, lj/cut/omp,
@ -48,7 +51,9 @@ explicitly in your input script, e.g. pair_style lj/cut/gpu. If the
suffix command is used with the appropriate style, you do not need to
modify your input script. The specified suffix (opt,omp,gpu,cuda) is
automatically appended whenever your input script command creates a
new "atom"_atom_style.html, "pair"_pair_style.html, "fix"_fix.html,
new "atom"_atom_style.html, "pair"_pair_style.html, "bond"_bond_style.html,
"angle"_angle_style.html, "dihedral"_dihedral_style.html,
"improper"_improper_style.html, "kspace"_kspace_style.html, "fix"_fix.html,
"compute"_compute.html, or "run"_run_style.html style. If the variant
version does not exist, the standard version is created.

8
src/.gitignore vendored
View File

@ -527,6 +527,14 @@ pair_lj_cut_coul_cut_gpu.cpp
pair_lj_cut_coul_cut_gpu.h
pair_lj_cut_coul_long_gpu.cpp
pair_lj_cut_coul_long_gpu.h
pair_lj_sdk.cpp
pair_lj_sdk.h
pair_lj_sdk_omp.cpp
pair_lj_sdk_omp.h
pair_lj_sdk_coul_long.cpp
pair_lj_sdk_coul_long.h
pair_lj_sdk_coul_long_omp.cpp
pair_lj_sdk_coul_long_omp.h
pair_lj_sf.cpp
pair_lj_sf.h
pair_lj_sf_omp.cpp

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

@ -15,6 +15,7 @@ if (test $1 = 1) then
cp pair_cdeam.cpp ..
cp pair_dipole_sf.cpp ..
cp pair_edip.cpp ..
cp pair_lj_sdk.cpp ..
cp pair_lj_sf.cpp ..
cp angle_cosine_shift.h ..
@ -27,9 +28,11 @@ if (test $1 = 1) then
cp fix_addtorque.h ..
cp fix_imd.h ..
cp fix_smd.h ..
cp lj_sdk_common.h ..
cp pair_cdeam.h ..
cp pair_dipole_sf.h ..
cp pair_edip.h ..
cp pair_lj_sdk.h ..
cp pair_lj_sf.h ..
elif (test $1 = 0) then
@ -47,6 +50,7 @@ elif (test $1 = 0) then
rm -f ../pair_cdeam.cpp
rm -f ../pair_dipole_sf.cpp
rm -f ../pair_edip.cpp
rm -f ../pair_lj_sdk.cpp
rm -f ../pair_lj_sf.cpp
rm -f ../angle_cosine_shift.h
@ -59,9 +63,11 @@ elif (test $1 = 0) then
rm -f ../fix_addtorque.h
rm -f ../fix_imd.h
rm -f ../fix_smd.h
rm -f ../lj_sdk_common.h
rm -f ../pair_cdeam.h
rm -f ../pair_dipole_sf.h
rm -f ../pair_edip.h
rm -f ../pair_lj_sdk.h
rm -f ../pair_lj_sf.h
fi

398
src/USER-MISC/lj_sdk_common.cpp
View File

@ -1,398 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#if 0
/* ----------------------------------------------------------------------
Common functionality for the SDK coarse grained MD potentials.
Contributing author: Axel Kohlmeyer <akohlmey@gmail.com>
------------------------------------------------------------------------- */
#include "lj_sdk_common.h"
#include "memory.h"
#include "stdlib.h"
#include "string.h"
#include "ctype.h"
#include "math.h"
using namespace LAMMPS_NS;
using namespace LJSDKParms;
#define SMALL 1.0e-6
/* ---------------------------------------------------------------------- */
PairLJSDKCommon::PairLJSDKCommon(class LAMMPS *lmp) : Pair(lmp)
{
respa_enable = 0;
single_enable = 0;
ftable = NULL;
cutsq = cut_lj = epsilon = sigma = lj1 = lj2 = lj3 = lj4 = offset = NULL;
allocated_coul = 0;
cut_lj_global = cut_coul_global = cut_coulsq_global = 0.0;
}
/* ---------------------------------------------------------------------- *
* clean up common arrays *
* ---------------------------------------------------------------------- */
PairLJSDKCommon::~PairLJSDKCommon() {
if (allocated) {
memory->destroy(setflag);
memory->destroy(cg_type);
memory->destroy(cutsq);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(offset);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
allocated = 0;
}
}
/* ---------------------------------------------------------------------- *
* allocate common arrays *
* ---------------------------------------------------------------------- */
void PairLJSDKCommon::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pairsdk:setflag");
memory->create(cg_type,n+1,n+1,"pairsdk:cg_type");
for (int i = 1; i <= n; i++) {
for (int j = i; j <= n; j++) {
setflag[i][j] = 0;
cg_type[i][j] = CG_NOT_SET;
}
}
memory->create(cutsq,n+1,n+1,"pairsdk:cutsq");
memory->create(epsilon,n+1,n+1,"pairsdk:epsilon");
memory->create(sigma,n+1,n+1,"pairsdk:sigma");
memory->create(offset,n+1,n+1,"pairsdk:offset");
memory->create(lj1,n+1,n+1,"pairsdk:lj1");
memory->create(lj2,n+1,n+1,"pairsdk:lj2");
memory->create(lj3,n+1,n+1,"pairsdk:lj3");
memory->create(lj4,n+1,n+1,"pairsdk:lj4");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
// arguments to the pair_style command (global version)
// args = cutoff (cutoff2)
void PairLJSDKCommon::settings(int narg, char **arg)
{
if ((narg < 1) || (narg > 3)) error->all(FLERR,"Illegal pair_style command");
cut_lj_global = force->numeric(arg[0]);
if (narg == 1) cut_coul_global = cut_lj_global;
else cut_coul_global = force->numeric(arg[1]);
// reset if we are in the no-charge version.
if (!allocated_coul) cut_coul_global=0.0;
cut_coulsq_global = cut_coul_global*cut_coul_global;
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++) {
for (j = i+1; j <= atom->ntypes; j++) {
if (setflag[i][j]) {
double cut = MAX(cut_lj_global,cut_coul_global);
cutsq[i][j] = cut*cut;
if (allocated_coul) {
cut_lj[i][j] = cut_lj_global;
}
}
}
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJSDKCommon::coeff(int narg, char **arg)
{
if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
int cg_type_one=find_cg_type(arg[2]);
if (cg_type_one == CG_NOT_SET) error->all(FLERR,"Error reading CG type flag.");
double epsilon_one = force->numeric(arg[3]);
double sigma_one = force->numeric(arg[4]);
double cut_lj_one = cut_lj_global;
double cut_coul_one = cut_coul_global;
if (narg >= 6) cut_lj_one = force->numeric(arg[5]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
cg_type[i][j] = cg_type_one;
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
setflag[i][j] = 1;
if (allocated_coul) {
cut_lj[i][j] = cut_lj_one;
cut_coul_one = MAX(cut_lj_one, cut_coul_global);
cutsq[i][j] = cut_coul_one * cut_coul_one;
} else {
cutsq[i][j] = cut_lj_one * cut_lj_one;
}
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJSDKCommon::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
error->all(FLERR,"for CG styles, epsilon and sigma need to be set explicitly for all pairs.");
}
const int cgt = cg_type[i][j];
if (cgt == CG_NOT_SET)
error->all(FLERR,"unrecognized LJ parameter flag");
lj1[i][j] = cg_prefact[cgt] * cg_pow1[cgt] * epsilon[i][j] * pow(sigma[i][j],cg_pow1[cgt]);
lj2[i][j] = cg_prefact[cgt] * cg_pow2[cgt] * epsilon[i][j] * pow(sigma[i][j],cg_pow2[cgt]);
lj3[i][j] = cg_prefact[cgt] * epsilon[i][j] * pow(sigma[i][j],cg_pow1[cgt]);
lj4[i][j] = cg_prefact[cgt] * epsilon[i][j] * pow(sigma[i][j],cg_pow2[cgt]);
double mycut = sqrt(cutsq[i][j]);
if (offset_flag) {
double ratio = sigma[i][j] / mycut;
offset[i][j] = cg_prefact[cgt] * epsilon[i][j] * (pow(ratio,cg_pow1[cgt]) - pow(ratio,cg_pow2[cgt]));
} else offset[i][j] = 0.0;
if (allocated_coul) {
mycut = MAX(cut_lj[i][j],cut_coul[i][j]);
cut[i][j] = mycut;
cut_ljsq[i][j]=cut_lj[i][j]*cut_lj[i][j];
cut_coulsq[i][j]=cut_coul[i][j]*cut_coul[i][j];
if (offset_flag) {
double ratio = sigma[i][j] / cut_lj[i][j];
offset[i][j] = cg_prefact[cgt] * epsilon[i][j] * (pow(ratio,cg_pow1[cgt]) - pow(ratio,cg_pow2[cgt]));
} else offset[i][j] = 0.0;
}
// make sure data is stored symmetrically
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
lj3[j][i] = lj3[i][j];
lj4[j][i] = lj4[i][j];
offset[j][i] = offset[i][j];
cg_type[j][i] = cg_type[i][j];
cut[j][i] = mycut;
if (allocated_coul) {
cut_lj[j][i]=cut_lj[i][j];
cut_ljsq[j][i]=cut_ljsq[i][j];
cut_coul[j][i]=cut_coul[i][j];
cut_coulsq[j][i]=cut_coulsq[i][j];
}
if (tail_flag) {
#if 1
error->all(FLERR,"tail correction not supported by CG potentials.");
#endif
}
return mycut;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJSDKCommon::write_restart(FILE *fp)
{
int i,j;
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&cg_type[i][j],sizeof(int),1,fp);
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
if (allocated_coul) {
fwrite(&cut_lj[i][j],sizeof(double),1,fp);
fwrite(&cut_coul[i][j],sizeof(double),1,fp);
}
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJSDKCommon::read_restart(FILE *fp)
{
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&cg_type[i][j],sizeof(int),1,fp);
fread(&epsilon[i][j],sizeof(double),1,fp);
fread(&sigma[i][j],sizeof(double),1,fp);
fread(&cut[i][j],sizeof(double),1,fp);
if(allocated_coul) {
fread(&cut_lj[i][j],sizeof(double),1,fp);
fread(&cut_coul[i][j],sizeof(double),1,fp);
}
}
MPI_Bcast(&cg_type[i][j],1,MPI_INT,0,world);
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
if (allocated_coul) {
MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world);
}
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJSDKCommon::write_restart_settings(FILE *fp)
{
fwrite(&cut_lj_global,sizeof(double),1,fp);
fwrite(&cut_coul_global,sizeof(double),1,fp);
fwrite(&kappa,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJSDKCommon::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
fread(&cut_lj_global,sizeof(double),1,fp);
fread(&cut_coul_global,sizeof(double),1,fp);
fread(&kappa,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&kappa,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
cut_coulsq_global = cut_coul_global*cut_coul_global;
}
/* ---------------------------------------------------------------------- */
double PairLJSDKCommon::memory_usage()
{
double bytes=Pair::memory_usage();
int n = atom->ntypes;
// setflag/cg_type
bytes += (n+1)*(n+1)*sizeof(int)*2;
// cut/cutsq/epsilon/sigma/offset/lj1/lj2/lj3/lj4
bytes += (n+1)*(n+1)*sizeof(double)*9;
return bytes;
}
/* ------------------------------------------------------------------------ */
double PairLJSDKCommon::eval_single(int coul_type, int i, int j, int itype, int jtype,
double rsq, double factor_coul, double factor_lj,
double &fforce)
{
double lj_force, lj_erg, coul_force, coul_erg;
lj_force=lj_erg=coul_force=coul_erg=0.0;
if (rsq < cut_ljsq[itype][jtype]) {
const int cgt = cg_type[itype][jtype];
const double cgpow1 = cg_pow1[cgt];
const double cgpow2 = cg_pow2[cgt];
const double cgpref = cg_prefact[cgt];
const double ratio = sigma[itype][jtype]/sqrt(rsq);
const double eps = epsilon[itype][jtype];
lj_force = cgpref*eps * (cgpow1*pow(ratio,cgpow1)
- cgpow2*pow(ratio,cgpow2))/rsq;
lj_erg = cgpref*eps * (pow(ratio,cgpow1) - pow(ratio,cgpow2));
}
if (rsq < cut_coul[itype][jtype]) {
if(coul_type == CG_COUL_LONG) {
error->all(FLERR,"single energy computation with long-range coulomb not supported by CG potentials.");
} else if ((coul_type == CG_COUL_CUT) || (coul_type == CG_COUL_DEBYE)) {
const double r2inv = 1.0/rsq;
const double rinv = sqrt(r2inv);
const double qscreen=exp(-kappa*sqrt(rsq));
coul_force = force->qqrd2e * atom->q[i]*atom->q[j]*rinv * qscreen * (kappa + rinv);
coul_erg = force->qqrd2e * atom->q[i]*atom->q[j]*rinv * qscreen;
// error->all(FLERR,"single energy computation with coulomb not supported by CG potentials.");
} else if (coul_type == CG_COUL_NONE) {
; // do nothing
} else {
error->all(FLERR,"unknown coulomb type with CG potentials.");
}
}
fforce = factor_lj*lj_force + factor_coul*coul_force;
return factor_lj*lj_erg + factor_coul*coul_erg;
}
#endif

89
src/USER-MISC/lj_sdk_common.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -10,11 +10,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#if 0
/* ----------------------------------------------------------------------
Common data for the Shinoda, DeVane, Klein (SDK) coars grain model
Contributing author: Axel Kohlmeyer (Temple U) <akohlmey@gmail.com>
Common data for the Shinoda, DeVane, Klein (SDK) coarse grain model
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#ifndef LMP_LJ_SDK_COMMON_H
@ -22,78 +21,24 @@
#include "string.h"
namespace LAMMPS_NS {
namespace LAMMPS_NS {
namespace LJSDKParms {
namespace LJSDKParms {
// LJ type flags. list of supported LJ exponent combinations
enum {LJ_NOT_SET=0, LJ9_6, LJ12_4, LJ12_6, NUM_LJ_TYPES};
// CG type flags. list of supported LJ exponent combinations
enum {CG_NOT_SET=0, CG_LJ9_6, CG_LJ12_4, CG_LJ12_6, NUM_CG_TYPES};
static int find_cg_type(const char *label) {
for (int i=0; i < NUM_CG_TYPES; ++i) {
if (strcmp(label,cg_type_list[i]) == 0) {
return i;
}
}
return CG_NOT_SET;
}
static const char * const cg_type_list[] = {"none", "lj9_6", "lj12_4", "lj12_6"};
static const double cg_prefact[] = {0.0, 6.75, 2.59807621135332, 4.0};
static const double cg_pow1[] = {0.0, 9.00, 12.0, 12.0};
static const double cg_pow2[] = {0.0, 6.00, 4.0, 6.0};
static int find_lj_type(const char *label,
const char * const * const list) {
for (int i=0; i < NUM_LJ_TYPES; ++i)
if (strcmp(label,list[i]) == 0) return i;
return LJ_NOT_SET;
}
class PairPairLJSDKCommon : public Pair {
public:
PairLJSDKCommon(class LAMMPS *);
virtual ~PairPairLJSDKCommon();
virtual void settings(int, char **);
virtual void coeff(int, char **);
virtual void init_style();
virtual void init_list(int, class NeighList *);
virtual double init_one(int, int);
virtual void write_restart(FILE *);
virtual void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual double memory_usage();
protected:
// coarse grain flags
int **cg_type;
// lennard jones parameters
double cut_lj_global, **cutsq, **cut_lj;
double **epsilon, **sigma;
double **lj1, **lj2, **lj3, **lj4, **offset;
// coulomb parameters
int allocated_coul; // 0/1 = whether coulomb arrays are allocated
double cut_coul_global, cut_coulsq_global, g_ewald;
// tables
double tabinnersq;
double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable;
double *etable,*detable,*ptable,*dptable,*vtable,*dvtable;
int ncoulshiftbits,ncoulmask;
// methods
virtual void allocate();
private:
// disable default constructor
PairLJSDKCommon();
}
static const char * const lj_type_list[] = {"none", "lj9_6", "lj12_4", "lj12_6"};
static const double lj_prefact[] = {0.0, 6.75, 2.59807621135332, 4.0};
static const double lj_pow1[] = {0.0, 9.00, 12.0, 12.0};
static const double lj_pow2[] = {0.0, 6.00, 4.0, 6.0};
}}
#endif
#endif

405
src/USER-MISC/pair_lj_sdk.cpp
View File

@ -11,22 +11,37 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#if 0
/* ----------------------------------------------------------------------
Shinoda, DeVane, Klein (SDK) model potential for coarse grained MD.
Plain version w/o charges.
Contributing author: Axel Kohlmeyer <akohlmey@gmail.com>
Contributing author: Axel Kohlmeyer (Temple U)
This style is a simplified re-implementation of the CG/CMM pair style
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_lj_sdk.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "lj_sdk_common.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace LJSDKParms;
/* ---------------------------------------------------------------------- */
PairLJSDK::PairLJSDK(LAMMPS *lmp) : PairCMMCommon(lmp)
PairLJSDK::PairLJSDK(LAMMPS *lmp) : Pair(lmp)
{
respa_enable = 0;
single_enable = 1;
@ -36,74 +51,376 @@ PairLJSDK::PairLJSDK(LAMMPS *lmp) : PairCMMCommon(lmp)
PairLJSDK::~PairLJSDK()
{
/* empty */ ;
if (allocated) {
memory->destroy(setflag);
memory->destroy(lj_type);
memory->destroy(cutsq);
memory->destroy(cut);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
memory->destroy(offset);
}
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- *
* the real compute work is done in the PairCMMCommon::eval_XXX<>() templates
* in the common PairCG class. Through using templates we can have one
* implementation for all CG varieties _and_ gain speed through having
* the compiler optimize away conditionals within the innerloops that
* can be predetermined outside the loop through instantiation of the
* different combination of template flags.
* ---------------------------------------------------------------------- */
void PairLJSDK::compute(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else {
evflag = vflag_fdotr = 0;
}
} else evflag = vflag_fdotr = 0;
if (evflag) {
if (eflag) {
if (force->newton_pair) {
return eval_verlet<1,1,1,CG_COUL_NONE>();
} else {
return eval_verlet<1,1,0,CG_COUL_NONE>();
}
if (force->newton_pair) eval<1,1,1>();
else eval<1,1,0>();
} else {
if (force->newton_pair) {
return eval_verlet<1,0,1,CG_COUL_NONE>();
} else {
return eval_verlet<1,0,0,CG_COUL_NONE>();
}
if (force->newton_pair) eval<1,0,1>();
else eval<1,0,0>();
}
} else {
if (force->newton_pair) {
return eval_verlet<0,0,1,CG_COUL_NONE>();
} else {
return eval_verlet<0,0,0,CG_COUL_NONE>();
if (force->newton_pair) eval<0,0,1>();
else eval<0,0,0>();
}
if (vflag_fdotr) virial_fdotr_compute();
}
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairLJSDK::eval()
{
int i,j,ii,jj,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r2inv,r6inv,forcelj,factor_lj;
evdwl = 0.0;
const double * const * const x = atom->x;
double * const * const f = atom->f;
const int * const type = atom->type;
const int nlocal = atom->nlocal;
const double * const special_lj = force->special_lj;
double fxtmp,fytmp,fztmp;
const int inum = list->inum;
const int * const ilist = list->ilist;
const int * const numneigh = list->numneigh;
const int * const * const firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
fxtmp=fytmp=fztmp=0.0;
const int itype = type[i];
const int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv;
fxtmp += delx*fpair;
fytmp += dely*fpair;
fztmp += delz*fpair;
if (NEWTON_PAIR || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (EFLAG) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype])
- offset[itype][jtype];
evdwl *= factor_lj;
}
if (EVFLAG) ev_tally(i,j,nlocal,NEWTON_PAIR,
evdwl,0.0,fpair,delx,dely,delz);
}
}
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJSDK::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(lj_type,n+1,n+1,"pair:lj_type");
for (int i = 1; i <= n; i++) {
for (int j = i; j <= n; j++) {
setflag[i][j] = 0;
lj_type[i][j] = LJ_NOT_SET;
}
}
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(lj1,n+1,n+1,"pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4");
memory->create(offset,n+1,n+1,"pair:offset");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJSDK::settings(int narg, char **arg)
{
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
cut_global = force->numeric(arg[0]);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i+1; j <= atom->ntypes; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
}
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJSDK::write_restart(FILE *fp)
void PairLJSDK::coeff(int narg, char **arg)
{
write_restart_settings(fp);
PairCMMCommon::write_restart(fp);
if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
int lj_type_one = find_lj_type(arg[2],lj_type_list);
if (lj_type_one == LJ_NOT_SET) error->all(FLERR,"Cannot parse LJ type flag.");
double epsilon_one = force->numeric(arg[2]);
double sigma_one = force->numeric(arg[3]);
double cut_one = cut_global;
if (narg == 5) cut_one = force->numeric(arg[4]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLJSDK::init_style()
{
neighbor->request(this);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJSDK::init_one(int i, int j)
{
if (setflag[i][j] == 0)
error->all(FLERR,"Mixing not supported for lj/sdk pair style");
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
if (offset_flag) {
error->all(FLERR,"Offset not supported for lj/sdk pair style");
double ratio = sigma[i][j] / cut[i][j];
offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
} else offset[i][j] = 0.0;
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
lj3[j][i] = lj3[i][j];
lj4[j][i] = lj4[i][j];
offset[j][i] = offset[i][j];
// compute I,J contribution to long-range tail correction
// count total # of atoms of type I and J via Allreduce
if (tail_flag) {
error->all(FLERR,"Tail flag not supported by lj/sdk pair style");
int *type = atom->type;
int nlocal = atom->nlocal;
double count[2],all[2];
count[0] = count[1] = 0.0;
for (int k = 0; k < nlocal; k++) {
if (type[k] == i) count[0] += 1.0;
if (type[k] == j) count[1] += 1.0;
}
MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
double sig2 = sigma[i][j]*sigma[i][j];
double sig6 = sig2*sig2*sig2;
double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
double rc6 = rc3*rc3;
double rc9 = rc3*rc6;
etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (sig6 - 3.0*rc6) / (9.0*rc9);
ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
}
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJSDK::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJSDK::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
PairCMMCommon::read_restart(fp);
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&epsilon[i][j],sizeof(double),1,fp);
fread(&sigma[i][j],sizeof(double),1,fp);
fread(&cut[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJSDK::write_restart_settings(FILE *fp)
{
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJSDK::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
fread(&cut_global,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
double PairLJSDK::single(int i, int j, int itype, int jtype, double rsq,
double factor_coul, double factor_lj, double &fforce)
double PairLJSDK::single(int, int, int itype, int jtype, double rsq,
double, double factor_lj, double &fforce)
{
return eval_single(CG_COUL_NONE,i,j,itype,jtype,rsq,factor_coul,factor_lj,fforce);
double r2inv,r6inv,forcelj,philj;
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
fforce = factor_lj*forcelj*r2inv;
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype];
return factor_lj*philj;
}
#endif
/* ---------------------------------------------------------------------- */
void *PairLJSDK::extract(char *str, int &dim)
{
dim = 2;
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
if (strcmp(str,"sigma") == 0) return (void *) sigma;
if (strcmp(str,"lj_type") == 0) return (void *) lj_type;
return NULL;
}

48
src/USER-MISC/pair_lj_sdk.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -11,8 +11,6 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#if 0
#ifdef PAIR_CLASS
PairStyle(lj/sdk,PairLJSDK)
@ -23,32 +21,44 @@ PairStyle(cg/cmm,PairLJSDK)
#ifndef LMP_PAIR_LJ_SDK_H
#define LMP_PAIR_LJ_SDK_H
#include "lj_sdk_common.h"
#include "pair.h"
namespace LAMMPS_NS {
class LAMMPS;
class PairLJSDK : public Pair {
public:
PairLJSDK(LAMMPS *);
virtual ~PairLJSDK();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
void *extract(char *, int &);
class PairLJSDK : PairLJSDKCommon {
protected:
int **lj_type; // type of lennard jones potential
public:
double **cut;
double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4,**offset;
PairLJSDK(class LAMMPS *);
virtual ~PairLJSDK();
double cut_global;
virtual void compute(int, int);
void allocate();
void write_restart(FILE *);
void read_restart(FILE *);
private:
template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void eval();
virtual double single(int, int, int, int, double, double, double, double &);
};
private:
// disable default constructor
PairLJSDK();
};
}
#endif
#endif
#endif