From 912d256f8e68bece12aacbd213536e25ed0eb1b5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 25 Oct 2011 23:02:07 -0400 Subject: [PATCH 1/3] improve documentation of suffix command --- doc/suffix.txt | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/doc/suffix.txt b/doc/suffix.txt index ec773f0eb7..5ee019a40d 100644 --- a/doc/suffix.txt +++ b/doc/suffix.txt @@ -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. From 997b9b04e571456c360c39cee522f1a733198bfc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 25 Oct 2011 23:04:03 -0400 Subject: [PATCH 2/3] add template for lj/sdk re-implementation --- src/USER-MISC/Install.sh | 6 + src/USER-MISC/lj_sdk_common.cpp | 398 ------------------------------- src/USER-MISC/lj_sdk_common.h | 89 ++----- src/USER-MISC/pair_lj_sdk.cpp | 405 ++++++++++++++++++++++++++++---- src/USER-MISC/pair_lj_sdk.h | 48 ++-- 5 files changed, 413 insertions(+), 533 deletions(-) delete mode 100644 src/USER-MISC/lj_sdk_common.cpp diff --git a/src/USER-MISC/Install.sh b/src/USER-MISC/Install.sh index d15156fdbd..6f17d9d806 100644 --- a/src/USER-MISC/Install.sh +++ b/src/USER-MISC/Install.sh @@ -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 diff --git a/src/USER-MISC/lj_sdk_common.cpp b/src/USER-MISC/lj_sdk_common.cpp deleted file mode 100644 index 60fbadeb0a..0000000000 --- a/src/USER-MISC/lj_sdk_common.cpp +++ /dev/null @@ -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 -------------------------------------------------------------------------- */ - -#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 diff --git a/src/USER-MISC/lj_sdk_common.h b/src/USER-MISC/lj_sdk_common.h index 899f86c5f3..9b151b8a35 100644 --- a/src/USER-MISC/lj_sdk_common.h +++ b/src/USER-MISC/lj_sdk_common.h @@ -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) + 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 diff --git a/src/USER-MISC/pair_lj_sdk.cpp b/src/USER-MISC/pair_lj_sdk.cpp index bc0e1e24fe..8cfdcf0089 100644 --- a/src/USER-MISC/pair_lj_sdk.cpp +++ b/src/USER-MISC/pair_lj_sdk.cpp @@ -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 + 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 +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; +} diff --git a/src/USER-MISC/pair_lj_sdk.h b/src/USER-MISC/pair_lj_sdk.h index c1c454468a..85d2da687a 100644 --- a/src/USER-MISC/pair_lj_sdk.h +++ b/src/USER-MISC/pair_lj_sdk.h @@ -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 void eval(); - virtual double single(int, int, int, int, double, double, double, double &); +}; - private: - // disable default constructor - PairLJSDK(); - }; } #endif #endif - - -#endif From 3c110b35da7cb887e9269b403ec5357b33a0e169 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 25 Oct 2011 23:06:16 -0400 Subject: [PATCH 3/3] more files to ignore --- src/.gitignore | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/.gitignore b/src/.gitignore index 67dd3f40ec..5a8ae74572 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -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