diff --git a/src/Makefile b/src/Makefile index c367e5dd55..88b7b09b7f 100755 --- a/src/Makefile +++ b/src/Makefile @@ -46,7 +46,7 @@ PACKAGE = asphere body class2 colloid coreshell dipole fld gpu granular kim \ python qeq reax replica rigid shock snap srd voronoi xtc PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars user-cuda \ - user-diffraction user-eff user-fep user-intel user-lb \ + user-diffraction user-drude user-eff user-fep user-intel user-lb \ user-misc user-molfile user-omp user-phonon user-qmmm user-qtb \ user-quip user-reaxc user-sph diff --git a/src/USER-DRUDE/README b/src/USER-DRUDE/README new file mode 100644 index 0000000000..e2d91a671e --- /dev/null +++ b/src/USER-DRUDE/README @@ -0,0 +1,21 @@ +This package implements Drude oscillator methods for simulating +polarizable systems in LAMMPS. The package provides the following +features: + +* tagging particles as cores or Drude particles +* thermostating the Drude oscillators at a distinct temperature +using Langevin or Nosé-Hoover thermostats +* computation of the atom and dipole temperatures +* damping induced dipole interactions using Thole's function + +See the file doc/drude_tutorial.html for getting started. + +There are auxiliary tools for using this package in tools/drude. + +There are example scripts for using this package in examples/USER/drude. + +The person who created this package is Alain Dequidt at the +Chemistry Institute of Clermont-Ferrand, Clermont University, France +(alain.dequidt at univ-bpclermont.fr). Contact him directly if you +have questions. Co-authors: Julien Devémy, Agilio Padua. + diff --git a/src/USER-DRUDE/compute_temp_drude.cpp b/src/USER-DRUDE/compute_temp_drude.cpp new file mode 100644 index 0000000000..97d264647d --- /dev/null +++ b/src/USER-DRUDE/compute_temp_drude.cpp @@ -0,0 +1,218 @@ +/* ---------------------------------------------------------------------- + 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 "mpi.h" +#include "stdlib.h" +#include "string.h" +#include "compute_temp_drude.h" +#include "atom.h" +#include "update.h" +#include "force.h" +#include "group.h" +#include "modify.h" +#include "fix.h" +#include "domain.h" +#include "lattice.h" +#include "memory.h" +#include "error.h" +#include "comm.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeTempDrude::ComputeTempDrude(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 3) error->all(FLERR,"Illegal compute temp command"); + + vector_flag = 1; + size_vector = 6; + extscalar = 0; + extvector = -1; + extlist = new int[6]; + extlist[0] = extlist[1] = 0; + extlist[2] = extlist[3] = extlist[4] = extlist[5] = 1; + tempflag = 0; // because does not compute a single temperature (scalar and vector) + + vector = new double[6]; + fix_drude = NULL; + id_temp = NULL; + temperature = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeTempDrude::~ComputeTempDrude() +{ + delete [] vector; + delete [] extlist; + delete [] id_temp; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempDrude::init() +{ + int ifix; + for (ifix = 0; ifix < modify->nfix; ifix++) + if (strcmp(modify->fix[ifix]->style,"drude") == 0) break; + if (ifix == modify->nfix) error->all(FLERR, "compute temp/drude requires fix drude"); + fix_drude = (FixDrude *) modify->fix[ifix]; + + if (!comm->ghost_velocity) + error->all(FLERR,"compute temp/drude requires ghost velocities. Use comm_modify vel yes"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempDrude::setup() +{ + dof_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempDrude::dof_compute() +{ + int nlocal = atom->nlocal; + int *type = atom->type; + int dim = domain->dimension; + int *drudetype = fix_drude->drudetype; + + fix_dof = 0; + for (int i = 0; i < modify->nfix; i++) + fix_dof += modify->fix[i]->dof(igroup); + + bigint dof_core_loc = 0, dof_drude_loc = 0; + for (int i = 0; i < nlocal; i++) { + if (atom->mask[i] & groupbit) { + if (drudetype[type[i]] == DRUDE_TYPE) // Non-polarizable atom + dof_drude_loc++; + else + dof_core_loc++; + } + } + dof_core_loc *= dim; + dof_drude_loc *= dim; + MPI_Allreduce(&dof_core_loc, &dof_core, 1, MPI_LMP_BIGINT, MPI_SUM, world); + MPI_Allreduce(&dof_drude_loc, &dof_drude, 1, MPI_LMP_BIGINT, MPI_SUM, world); + dof_core -= fix_dof; + vector[2] = dof_core; + vector[3] = dof_drude; +} + +/* ---------------------------------------------------------------------- */ + +int ComputeTempDrude::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"temp") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + delete [] id_temp; + int n = strlen(arg[1]) + 1; + id_temp = new char[n]; + strcpy(id_temp,arg[1]); + + int icompute = modify->find_compute(id_temp); + if (icompute < 0) + error->all(FLERR,"Could not find fix_modify temperature ID"); + temperature = modify->compute[icompute]; + + if (temperature->tempflag == 0) + error->all(FLERR, + "Fix_modify temperature ID does not compute temperature"); + if (temperature->igroup != igroup && comm->me == 0) + error->warning(FLERR,"Group for fix_modify temp != fix group"); + return 2; + } + return 0; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempDrude::compute_vector() +{ + invoked_vector = update->ntimestep; + + int nlocal = atom->nlocal; + int *mask = atom->mask; + int *type = atom->type; + double *rmass = atom->rmass, *mass = atom->mass; + double **v = atom->v; + tagint *drudeid = fix_drude->drudeid; + int *drudetype = fix_drude->drudetype; + int dim = domain->dimension; + double mvv2e = force->mvv2e, kb = force->boltz; + + double mcore, mdrude; + double ecore, edrude; + double *vcore, *vdrude; + double kineng_core_loc = 0., kineng_drude_loc = 0.; + for (int i=0; iremove_bias(i, vcore); + for (int k=0; krestore_bias(i, vcore); + if (rmass) mcore = rmass[i]; + else mcore = mass[type[i]]; + kineng_core_loc += mcore * ecore; + } else { // CORE_TYPE + int j = atom->map(drudeid[i]); + if (rmass) { + mcore = rmass[i]; + mdrude = rmass[j]; + } else { + mcore = mass[type[i]]; + mdrude = mass[type[j]]; + } + double mtot_inv = 1. / (mcore + mdrude); + ecore = 0.; + edrude = 0.; + vcore = v[i]; + vdrude = v[j]; + if (temperature) { + temperature->remove_bias(i, vcore); + temperature->remove_bias(j, vdrude); + } + for (int k=0; krestore_bias(i, vcore); + temperature->restore_bias(j, vdrude); + } + kineng_core_loc += mtot_inv * ecore; + kineng_drude_loc += mtot_inv * mcore * mdrude * edrude; + } + } + } + + if (dynamic) dof_compute(); + kineng_core_loc *= 0.5 * mvv2e; + kineng_drude_loc *= 0.5 * mvv2e; + MPI_Allreduce(&kineng_core_loc,&kineng_core,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&kineng_drude_loc,&kineng_drude,1,MPI_DOUBLE,MPI_SUM,world); + temp_core = 2.0 * kineng_core / (dof_core * kb); + temp_drude = 2.0 * kineng_drude / (dof_drude * kb); + vector[0] = temp_core; + vector[1] = temp_drude; + vector[4] = kineng_core; + vector[5] = kineng_drude; +} + diff --git a/src/USER-DRUDE/compute_temp_drude.h b/src/USER-DRUDE/compute_temp_drude.h new file mode 100644 index 0000000000..834d1d7519 --- /dev/null +++ b/src/USER-DRUDE/compute_temp_drude.h @@ -0,0 +1,63 @@ +/* -*- 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(temp/drude,ComputeTempDrude) + +#else + +#ifndef LMP_COMPUTE_TEMP_DRUDE_H +#define LMP_COMPUTE_TEMP_DRUDE_H + +#include "compute.h" +#include "fix_drude.h" + +namespace LAMMPS_NS { + +class ComputeTempDrude : public Compute { + public: + ComputeTempDrude(class LAMMPS *, int, char **); + ~ComputeTempDrude(); + void init(); + void setup(); + void compute_vector(); + int modify_param(int, char **); + + private: + int fix_dof; + FixDrude * fix_drude; + char *id_temp; + class Compute *temperature; + bigint dof_core, dof_drude; + double kineng_core, kineng_drude; + double temp_core, temp_drude; + + void dof_compute(); + +}; + +} + +#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. + +*/ diff --git a/src/USER-DRUDE/fix_drude.cpp b/src/USER-DRUDE/fix_drude.cpp new file mode 100644 index 0000000000..9b95c15060 --- /dev/null +++ b/src/USER-DRUDE/fix_drude.cpp @@ -0,0 +1,496 @@ +/* ---------------------------------------------------------------------- + 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 "stdlib.h" +#include "fix_drude.h" +#include "atom.h" +#include "comm.h" +#include "modify.h" +#include "error.h" +#include "memory.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +FixDrude *FixDrude::sptr = NULL; + +/* ---------------------------------------------------------------------- */ + +FixDrude::FixDrude(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg != 3 + atom->ntypes) error->all(FLERR,"Illegal fix drude command"); + + comm_border = 1; // drudeid + special_alter_flag = 1; + create_attribute = 1; + + memory->create(drudetype, atom->ntypes+1, "fix_drude::drudetype"); + for (int i=3; iall(FLERR, "Illegal fix drude command"); + } + + drudeid = NULL; + grow_arrays(atom->nmax); + atom->add_callback(0); + atom->add_callback(1); + atom->add_callback(2); + + // one-time assignment of Drude partners + + build_drudeid(); + + // set rebuildflag = 0 to indicate special lists have never been rebuilt + + rebuildflag = 0; +} + +/* ---------------------------------------------------------------------- */ + +FixDrude::~FixDrude() +{ + atom->delete_callback(id,2); + atom->delete_callback(id,1); + atom->delete_callback(id,0); + memory->destroy(drudetype); + memory->destroy(drudeid); +} + +/* ---------------------------------------------------------------------- */ + +void FixDrude::init() +{ + int count = 0; + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"drude") == 0) count++; + if (count > 1) error->all(FLERR,"More than one fix drude"); + + if (!rebuildflag) rebuild_special(); +} + +/* ---------------------------------------------------------------------- */ + +int FixDrude::setmask() +{ + int mask = 0; + return mask; +} + +/* ---------------------------------------------------------------------- + look in bond lists for Drude partner tags and fill drudeid +------------------------------------------------------------------------- */ +void FixDrude::build_drudeid(){ + int nlocal = atom->nlocal; + int *type = atom->type; + + std::vector drude_vec; // list of my Drudes' tags + std::vector core_drude_vec; + partner_set = new std::set[nlocal]; // Temporary sets of bond partner tags + + sptr = this; + // Build list of my atoms' bond partners + for (int i=0; inum_bond[i]; k++){ + core_drude_vec.push_back(atom->tag[i]); + core_drude_vec.push_back(atom->bond_atom[i][k]); + } + } + // Loop on procs to fill my atoms' sets of bond partners + comm->ring(core_drude_vec.size(), sizeof(tagint), + (char *) core_drude_vec.data(), + 4, ring_build_partner, NULL, 1); + + // Build the list of my Drudes' tags + // The only bond partners of a Drude particle is its core, + // so fill drudeid for my Drudes. + for (int i=0; itag[i]); + drudeid[i] = *partner_set[i].begin(); // only one 1-2 neighbor, the core + } + } + // At this point each of my Drudes knows its core. + // Send my list of Drudes to other procs and myself + // so that each core finds its Drude. + comm->ring(drude_vec.size(), sizeof(tagint), + (char *) drude_vec.data(), + 3, ring_search_drudeid, NULL, 1); + delete [] partner_set; +} + +/* ---------------------------------------------------------------------- + * when receive buffer, build the set of received Drude tags. + * Look in my cores' bond partner tags if there is a Drude tag. + * If so fill this core's dureid. +------------------------------------------------------------------------- */ +void FixDrude::ring_search_drudeid(int size, char *cbuf){ + // Search for the drude partner of my cores + Atom *atom = sptr->atom; + int nlocal = atom->nlocal; + int *type = atom->type; + std::set *partner_set = sptr->partner_set; + tagint *drudeid = sptr->drudeid; + int *drudetype = sptr->drudetype; + + tagint *first = (tagint *) cbuf; + tagint *last = first + size; + std::set drude_set(first, last); + std::set::iterator it; + + for (int i=0; i 0) continue; + for (it = partner_set[i].begin(); it != partner_set[i].end(); it++) { // Drude-core are 1-2 neighbors + if (drude_set.count(*it) > 0){ + drudeid[i] = *it; + break; + } + } + } +} + +/* ---------------------------------------------------------------------- + * buffer contains bond partners. Look for my atoms and add their partner's + * tag in its set of bond partners. +------------------------------------------------------------------------- */ +void FixDrude::ring_build_partner(int size, char *cbuf){ + // Add partners from incoming list + Atom *atom = sptr->atom; + int nlocal = atom->nlocal; + std::set *partner_set = sptr->partner_set; + tagint *it = (tagint *) cbuf; + tagint *last = it + size; + + while (it < last) { + int j = atom->map(*it); + if (j >= 0 && j < nlocal) + partner_set[j].insert(*(it+1)); + j = atom->map(*(it+1)); + if (j >= 0 && j < nlocal) + partner_set[j].insert(*it); + it += 2; + } +} + + +/* ---------------------------------------------------------------------- + allocate atom-based array for drudeid +------------------------------------------------------------------------- */ + +void FixDrude::grow_arrays(int nmax) +{ + memory->grow(drudeid,nmax,"fix_drude:drudeid"); +} + +/* ---------------------------------------------------------------------- + copy values within local atom-based array +------------------------------------------------------------------------- */ + +void FixDrude::copy_arrays(int i, int j, int delflag) +{ + drudeid[j] = drudeid[i]; +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based array for exchange with another proc +------------------------------------------------------------------------- */ + +int FixDrude::pack_exchange(int i, double *buf) +{ + int m = 0; + buf[m++] = ubuf(drudeid[i]).d; + return m; +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based array from exchange with another proc +------------------------------------------------------------------------- */ + +int FixDrude::unpack_exchange(int nlocal, double *buf) +{ + int m = 0; + drudeid[nlocal] = (tagint) ubuf(buf[m++]).i; + return m; +} + +/* ---------------------------------------------------------------------- + pack values for border communication at re-neighboring +------------------------------------------------------------------------- */ + +int FixDrude::pack_border(int n, int *list, double *buf) +{ + int m = 0; + for (int i=0; inlocal; + int **nspecial = atom->nspecial; + tagint **special = atom->special; + int *type = atom->type; + + // Make sure that drude partners know each other + //build_drudeid(); + + // Log info + if (comm->me == 0) { + if (screen) fprintf(screen, "Rebuild special list taking Drude particles into account\n"); + if (logfile) fprintf(logfile, "Rebuild special list taking Drude particles into account\n"); + } + int nspecmax, nspecmax_old, nspecmax_loc; + nspecmax_loc = 0; + for (int i=0; ime == 0) { + if (screen) fprintf(screen, "Old max number of 1-2 to 1-4 neighbors: %d\n", nspecmax_old); + if (logfile) fprintf(logfile, "Old max number of 1-2 to 1-4 neighbors: %d\n", nspecmax_old); + } + + // Build lists of drude and core-drude pairs + std::vector drude_vec, core_drude_vec, core_special_vec; + for (int i=0; itag[i]); + } else if (drudetype[type[i]] == CORE_TYPE){ + core_drude_vec.push_back(atom->tag[i]); + core_drude_vec.push_back(drudeid[i]); + } + } + // Remove Drude particles from the special lists of each proc + comm->ring(drude_vec.size(), sizeof(tagint), + (char *) drude_vec.data(), + 9, ring_remove_drude, NULL, 1); + // Add back Drude particles in the lists just after their core + comm->ring(core_drude_vec.size(), sizeof(tagint), + (char *) core_drude_vec.data(), + 10, ring_add_drude, NULL, 1); + + // Check size of special list + nspecmax_loc = 0; + for (int i=0; ime == 0) { + if (screen) fprintf(screen, "New max number of 1-2 to 1-4 neighbors: %d (+%d)\n", nspecmax, nspecmax - nspecmax_old); + if (logfile) fprintf(logfile, "New max number of 1-2 to 1-4 neighbors: %d (+%d)\n", nspecmax, nspecmax - nspecmax_old); + } + if (atom->maxspecial < nspecmax) { + char str[1024]; + sprintf(str, "Not enough space in special: special_bonds extra should be at least %d", nspecmax - nspecmax_old); + error->all(FLERR, str); + } + + // Build list of cores' special lists to communicate to ghost drude particles + for (int i=0; itag[i]); + core_special_vec.push_back((tagint) nspecial[i][0]); + core_special_vec.push_back((tagint) nspecial[i][1]); + core_special_vec.push_back((tagint) nspecial[i][2]); + for (int j=1; jring(core_special_vec.size(), sizeof(tagint), + (char *) core_special_vec.data(), + 11, ring_copy_drude, NULL, 1); +} + +/* ---------------------------------------------------------------------- + * When receive buffer, build a set of drude tags, look into my atoms' + * special list if some tags are drude particles. If so, remove it. +------------------------------------------------------------------------- */ +void FixDrude::ring_remove_drude(int size, char *cbuf){ + // Remove all drude particles from special list + Atom *atom = sptr->atom; + int nlocal = atom->nlocal; + int **nspecial = atom->nspecial; + tagint **special = atom->special; + int *type = atom->type; + tagint *first = (tagint *) cbuf; + tagint *last = first + size; + std::set drude_set(first, last); + int *drudetype = sptr->drudetype; + + for (int i=0; i 0) { // I identify a drude in my special list, remove it + // left shift + nspecial[i][2]--; + for (int k=j; k drude tag. + * Loop on my atoms' special list to find core tags. Insert their Drude + * particle if they have one. +------------------------------------------------------------------------- */ +void FixDrude::ring_add_drude(int size, char *cbuf){ + // Assume special array size is big enough + // Add all particle just after their core in the special list + Atom *atom = sptr->atom; + int nlocal = atom->nlocal; + int **nspecial = atom->nspecial; + tagint **special = atom->special; + int *type = atom->type; + tagint *drudeid = sptr->drudeid; + int *drudetype = sptr->drudetype; + + tagint *first = (tagint *) cbuf; + tagint *last = first + size; + std::map core_drude_map; + + tagint *it = first; + while (it < last) { + tagint core_tag = *it; + it++; + core_drude_map[core_tag] = *it; + it++; + } + + for (int i=0; itag[i]) > 0) { // I identify myself as a core, add my own drude + // right shift + for (int k=nspecial[i][2]-1; k>=0; k--) + special[i][k+1] = special[i][k]; + special[i][0] = drudeid[i]; + nspecial[i][0]++; + nspecial[i][1]++; + nspecial[i][2]++; + } + for (int j=0; j 0) { // I identify a core in my special list, add his drude + // right shift + for (int k=nspecial[i][2]-1; k>j; k--) + special[i][k+1] = special[i][k]; + special[i][j+1] = core_drude_map[special[i][j]]; + nspecial[i][2]++; + if (j < nspecial[i][1]) { + nspecial[i][1]++; + if (j < nspecial[i][0]) nspecial[i][0]++; + } + j++; + } + } + } +} + +/* ---------------------------------------------------------------------- + * When receive buffer, build a map core tag -> pointer to special info + * in the buffer. Loop on my Drude particles and copy their special + * info from that of their core if the latter is found in the map. +------------------------------------------------------------------------- */ +void FixDrude::ring_copy_drude(int size, char *cbuf){ + // Copy special list of drude from its core (except itself) + Atom *atom = sptr->atom; + int nlocal = atom->nlocal; + int **nspecial = atom->nspecial; + tagint **special = atom->special; + int *type = atom->type; + tagint *drudeid = sptr->drudeid; + int *drudetype = sptr->drudetype; + + tagint *first = (tagint *) cbuf; + tagint *last = first + size; + std::map core_special_map; + + tagint *it = first; + while (it < last) { + tagint core_tag = *it; + it++; + core_special_map[core_tag] = it; + it += 2; + it += (int) *it; + } + + for (int i=0; i 0) { // My core is in this list, copy its special info + it = core_special_map[drudeid[i]]; + nspecial[i][0] = (int) *it; + it++; + nspecial[i][1] = (int) *it; + it++; + nspecial[i][2] = (int) *it; + it++; + special[i][0] = drudeid[i]; + for (int k=1; ktype[i]] != NOPOL_TYPE){ + if (atom->nspecial[i] ==0) error->all(FLERR, "Polarizable atoms cannot be inserted with special lists info from the molecule template"); + drudeid[i] = atom->special[i][0]; // Drude partner should be at first place in the special list + } else { + drudeid[i] = 0; + } +} + diff --git a/src/USER-DRUDE/fix_drude.h b/src/USER-DRUDE/fix_drude.h new file mode 100644 index 0000000000..ca2fc9fdb9 --- /dev/null +++ b/src/USER-DRUDE/fix_drude.h @@ -0,0 +1,69 @@ +/* -*- 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 FIX_CLASS + +FixStyle(drude,FixDrude) + +#else + +#ifndef LMP_FIX_DRUDE_H +#define LMP_FIX_DRUDE_H + +#include "fix.h" +#include + +#define NOPOL_TYPE 0 +#define CORE_TYPE 1 +#define DRUDE_TYPE 2 + +namespace LAMMPS_NS { + +class FixDrude : public Fix { + public: + int * drudetype; + tagint * drudeid; + bool is_reduced; + + FixDrude(class LAMMPS *, int, char **); + virtual ~FixDrude(); + int setmask(); + void init(); + + void grow_arrays(int nmax); + void copy_arrays(int i, int j, int delflag); + void set_arrays(int i); + int pack_exchange(int i, double *buf); + int unpack_exchange(int nlocal, double *buf); + int pack_border(int n, int *list, double *buf); + int unpack_border(int n, int first, double *buf); + +private: + int rebuildflag; + static FixDrude *sptr; + std::set * partner_set; + + void build_drudeid(); + static void ring_search_drudeid(int size, char *cbuf); + static void ring_build_partner(int size, char *cbuf); + void rebuild_special(); + static void ring_remove_drude(int size, char *cbuf); + static void ring_add_drude(int size, char *cbuf); + static void ring_copy_drude(int size, char *cbuf); +}; + +} + +#endif +#endif + diff --git a/src/USER-DRUDE/fix_drude_transform.cpp b/src/USER-DRUDE/fix_drude_transform.cpp new file mode 100644 index 0000000000..6d00b35695 --- /dev/null +++ b/src/USER-DRUDE/fix_drude_transform.cpp @@ -0,0 +1,317 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/** Fix Drude Transform ******************************************************/ +#include "math.h" +#include "fix_drude_transform.h" +#include "atom.h" +#include "domain.h" +#include "comm.h" +#include "error.h" +#include "modify.h" +#include "force.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ +template +FixDrudeTransform::FixDrudeTransform(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), mcoeff(NULL) +{ + if (narg != 3) error->all(FLERR,"Illegal fix drude/transform command"); + comm_forward = 9; + fix_drude = NULL; +} + +/* ---------------------------------------------------------------------- */ +template +FixDrudeTransform::~FixDrudeTransform() +{ + if (mcoeff) delete [] mcoeff; +} + +/* ---------------------------------------------------------------------- */ +template +void FixDrudeTransform::init() +{ + int ifix; + for (ifix = 0; ifix < modify->nfix; ifix++) + if (strcmp(modify->fix[ifix]->style,"drude") == 0) break; + if (ifix == modify->nfix) error->all(FLERR, "fix drude/transform requires fix drude"); + fix_drude = (FixDrude *) modify->fix[ifix]; +} + +/* ---------------------------------------------------------------------- */ +template +int FixDrudeTransform::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + return mask; +} + +/* ---------------------------------------------------------------------- */ +template +void FixDrudeTransform::setup(int) { + int nlocal = atom->nlocal; + int ntypes = atom->ntypes; + int * type = atom->type; + double * rmass = atom->rmass, * mass = atom->mass; + tagint * drudeid = fix_drude->drudeid; + int * drudetype = fix_drude->drudetype; + + if (!rmass) { + if (!mcoeff) mcoeff = new double[ntypes+1]; + double mcoeff_loc[ntypes+1]; + for (int itype=0; itype<=ntypes; itype++) mcoeff_loc[itype] = 2.; // an impossible value: mcoeff is at most 1. + for (int i=0; imap(drudeid[i]); + // i is drude, j is core + if (mcoeff_loc[type[i]] < 1.5) { // already done + if (mcoeff_loc[type[j]] > 1.5){ // not yet done ?? + error->all(FLERR,"There must be one Drude type per core type");} + continue; + } + mcoeff_loc[type[i]] = mass[type[i]] / (mass[type[i]] + mass[type[j]]); + mcoeff_loc[type[j]] = -mass[type[i]] / mass[type[j]]; + } + } + + MPI_Allreduce(mcoeff_loc, mcoeff, ntypes+1, MPI_DOUBLE, MPI_MIN, world); + // mcoeff is 2 for non polarizable + // 0 < mcoeff < 1 for drude + // mcoeff < 0 for core + } +} + +/* ---------------------------------------------------------------------- */ +namespace LAMMPS_NS { // required for specialization +template <> +void FixDrudeTransform::initial_integrate(int){ + comm->forward_comm_fix(this); + real_to_reduced(); + //comm->forward_comm_fix(this); // Normally not needed +} + +template <> +void FixDrudeTransform::final_integrate(){ + comm->forward_comm_fix(this); + real_to_reduced(); + //comm->forward_comm_fix(this); // Normally not needed +} + +template <> +void FixDrudeTransform::initial_integrate(int){ + comm->forward_comm_fix(this); + reduced_to_real(); + //comm->forward_comm_fix(this); // Normally not needed +} + +template <> +void FixDrudeTransform::final_integrate(){ + comm->forward_comm_fix(this); + reduced_to_real(); + //comm->forward_comm_fix(this); // Normally not needed +} + +} // end of namespace + +/* ---------------------------------------------------------------------- */ +template +void FixDrudeTransform::real_to_reduced() +{ + int nlocal = atom->nlocal; + int ntypes = atom->ntypes; + int dim = domain->dimension; + int * mask = atom->mask, * type = atom->type; + double ** x = atom->x, ** v = atom->v, ** f = atom->f; + double * rmass = atom->rmass, * mass = atom->mass; + double mcore, mdrude, coeff; + int icore, idrude; + tagint * drudeid = fix_drude->drudeid; + int * drudetype = fix_drude->drudetype; + + if (!rmass) { // TODO: maybe drudetype can be used instead? + for (int itype=1; itype<=ntypes; itype++) + if (mcoeff[itype] < 1.5) mass[itype] *= 1. - mcoeff[itype]; + } + for (int i=0; iclosest_image(i, atom->map(drudeid[i])); + } + } + for (int i=0; iis_reduced = true; +} + +/* ---------------------------------------------------------------------- */ +template +void FixDrudeTransform::reduced_to_real() +{ + int nlocal = atom->nlocal; + int ntypes = atom->ntypes; + int dim = domain->dimension; + int * mask = atom->mask, * type = atom->type; + double ** x = atom->x, ** v = atom->v, ** f = atom->f; + double * rmass = atom->rmass, * mass = atom->mass; + double mcore, mdrude, coeff; + int icore, idrude; + tagint * drudeid = fix_drude->drudeid; + int * drudetype = fix_drude->drudetype; + + for (int i=0; i 1.5 ? + double s = sqrt(1. - mass[type[idrude]] / mass[type[icore]]); + mass[type[idrude]] = 0.5 * mass[type[icore]] * (1. - s); + mdrude = mass[type[idrude]]; + mass[type[icore]] -= mdrude; + mcore = mass[type[icore]]; + mcoeff[type[icore]] = mdrude / (mcore + mdrude); + } + coeff = mcoeff[type[idrude]]; + } + for (int k=0; ktag[(int) drudeid[i]]; + } + } + if (!rmass) { + for (int itype=1; itype<=ntypes; itype++) + if (mcoeff[itype] < 1.5) mass[itype] /= 1. - mcoeff[itype]; + } + fix_drude->is_reduced = false; +} + +/* ---------------------------------------------------------------------- */ +template +int FixDrudeTransform::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) +{ + double ** x = atom->x, ** v = atom->v, ** f = atom->f; + int * type = atom->type, * drudetype = fix_drude->drudetype; + double dx,dy,dz; + int dim = domain->dimension; + int m = 0; + for (int i=0; iis_reduced && drudetype[type[j]] == DRUDE_TYPE)) { + for (int k=0; ktriclinic != 0) { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy; + dy = pbc[1]*domain->yprd; + if (dim == 3) { + dx += + pbc[4]*domain->xz; + dy += pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + } + else { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + if (dim == 3) + dz = pbc[2]*domain->zprd; + } + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + if (dim == 3) + buf[m++] = x[j][2] + dz; + } + for (int k=0; k +void FixDrudeTransform::unpack_forward_comm(int n, int first, double *buf) +{ + double ** x = atom->x, ** v = atom->v, ** f = atom->f; + int dim = domain->dimension; + int m = 0; + int last = first + n; + for (int i=first; i; +template class FixDrudeTransform; + diff --git a/src/USER-DRUDE/fix_drude_transform.h b/src/USER-DRUDE/fix_drude_transform.h new file mode 100644 index 0000000000..41f4e147ed --- /dev/null +++ b/src/USER-DRUDE/fix_drude_transform.h @@ -0,0 +1,52 @@ +/* -*- 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 FIX_CLASS + +FixStyle(drude/transform/direct,FixDrudeTransform) +FixStyle(drude/transform/inverse,FixDrudeTransform) + +#else + +#ifndef LMP_FIX_DRUDE_TRANSFORM_H +#define LMP_FIX_DRUDE_TRANSFORM_H + +#include "fix.h" +#include "fix_drude.h" + +namespace LAMMPS_NS { + +template +class FixDrudeTransform : public Fix { + public: + FixDrudeTransform(class LAMMPS *, int, char **); + ~FixDrudeTransform(); + int setmask(); + void init(); + void setup(int vflag); + void reduced_to_real(); + void real_to_reduced(); + void initial_integrate(int vflag); + void final_integrate(); + int pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc); + void unpack_forward_comm(int n, int first, double *buf); + protected: + double * mcoeff; + FixDrude * fix_drude; +}; + +} + +#endif +#endif + diff --git a/src/USER-DRUDE/fix_langevin_drude.cpp b/src/USER-DRUDE/fix_langevin_drude.cpp new file mode 100644 index 0000000000..041e7b79ae --- /dev/null +++ b/src/USER-DRUDE/fix_langevin_drude.cpp @@ -0,0 +1,420 @@ +/* ---------------------------------------------------------------------- + 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 "stdlib.h" +#include "math.h" +#include "fix_langevin_drude.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "input.h" +#include "variable.h" +#include "random_mars.h" +#include "group.h" +#include "update.h" +#include "modify.h" +#include "compute.h" +#include "error.h" +#include "domain.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{NOBIAS,BIAS}; +enum{CONSTANT,EQUAL}; + +/* ---------------------------------------------------------------------- */ + +FixLangevinDrude::FixLangevinDrude(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 9) error->all(FLERR,"Illegal fix langevin/drude command"); + // TODO add option for tally + + // Langevin thermostat should be applied every step + nevery = 1; + vector_flag = 1; + global_freq = nevery; + extvector = 0; + size_vector = 6; + comm_reverse = 3; + //extscalar = 1; + + // core temperature + tstr_core = NULL; + if (strstr(arg[3],"v_") == arg[3]) { + int n = strlen(&arg[3][2]) + 1; + tstr_core = new char[n]; + strcpy(tstr_core,&arg[3][2]); + tstyle_core = EQUAL; + } else { + t_start_core = force->numeric(FLERR,arg[3]); + t_target_core = t_start_core; + tstyle_core = CONSTANT; + } + t_period_core = force->numeric(FLERR,arg[4]); + int seed_core = force->inumeric(FLERR,arg[5]); + + // drude temperature + tstr_drude = NULL; + if (strstr(arg[7],"v_") == arg[6]) { + int n = strlen(&arg[6][2]) + 1; + tstr_drude = new char[n]; + strcpy(tstr_drude,&arg[6][2]); + tstyle_drude = EQUAL; + } else { + t_start_drude = force->numeric(FLERR,arg[6]); + t_target_drude = t_start_drude; + tstyle_drude = CONSTANT; + } + t_period_drude = force->numeric(FLERR,arg[7]); + int seed_drude = force->inumeric(FLERR,arg[8]); + + // error checks + if (t_period_core <= 0.0) + error->all(FLERR,"Fix langevin/drude period must be > 0.0"); + if (seed_core <= 0) error->all(FLERR,"Illegal langevin/drude seed"); + if (t_period_drude <= 0.0) + error->all(FLERR,"Fix langevin/drude period must be > 0.0"); + if (seed_drude <= 0) error->all(FLERR,"Illegal langevin/drude seed"); + + random_core = new RanMars(lmp,seed_core); + random_drude = new RanMars(lmp,seed_drude); + + int iarg = 9; + zero = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"zero") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin/drude command"); + if (strcmp(arg[iarg+1],"no") == 0) zero = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) zero = 1; + else error->all(FLERR,"Illegal fix langevin/drude command"); + iarg += 2; + } else error->all(FLERR,"Illegal fix langevin/drude command"); + } + + tflag = 0; // no external compute/temp is specified yet (for bias) + energy = 0.; + fix_drude = NULL; + temperature = NULL; + id_temp = NULL; +} + +/* ---------------------------------------------------------------------- */ + +FixLangevinDrude::~FixLangevinDrude() +{ + delete random_core; + delete [] tstr_core; + delete random_drude; + delete [] tstr_drude; +} + +/* ---------------------------------------------------------------------- */ + +int FixLangevinDrude::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixLangevinDrude::init() +{ + // check variable-style target core temperature + if (tstr_core) { + tvar_core = input->variable->find(tstr_core); + if (tvar_core < 0) + error->all(FLERR,"Variable name for fix langevin/drude does not exist"); + if (input->variable->equalstyle(tvar_core)) tstyle_core = EQUAL; + else error->all(FLERR,"Variable for fix langevin/drude is invalid style"); + } + + // check variable-style target drude temperature + if (tstr_drude) { + tvar_drude = input->variable->find(tstr_drude); + if (tvar_drude < 0) + error->all(FLERR,"Variable name for fix langevin/drude does not exist"); + if (input->variable->equalstyle(tvar_drude)) tstyle_drude = EQUAL; + else error->all(FLERR,"Variable for fix langevin/drude is invalid style"); + } + + int ifix; + for (ifix = 0; ifix < modify->nfix; ifix++) + if (strcmp(modify->fix[ifix]->style,"drude") == 0) break; + if (ifix == modify->nfix) error->all(FLERR, "fix langevin/drude requires fix drude"); + fix_drude = (FixDrude *) modify->fix[ifix]; +} + +/* ---------------------------------------------------------------------- */ + +void FixLangevinDrude::setup(int vflag) +{ + if (!strstr(update->integrate_style,"verlet")) + error->all(FLERR,"RESPA style not compatible with fix langevin/drude"); + if (!comm->ghost_velocity) + error->all(FLERR,"fix langevin/drude requires ghost velocities. Use comm_modify vel yes"); + + if (zero) { + int *mask = atom->mask; + int nlocal = atom->nlocal; + int *drudetype = fix_drude->drudetype; + int *type = atom->type; + bigint ncore_loc = 0; + for (int i=0; iall(FLERR,"Illegal fix_modify command"); + delete [] id_temp; + int n = strlen(arg[1]) + 1; + id_temp = new char[n]; + strcpy(id_temp,arg[1]); + + int icompute = modify->find_compute(id_temp); + if (icompute < 0) + error->all(FLERR,"Could not find fix_modify temperature ID"); + temperature = modify->compute[icompute]; + + if (temperature->tempflag == 0) + error->all(FLERR, + "Fix_modify temperature ID does not compute temperature"); + if (temperature->igroup != igroup && comm->me == 0) + error->warning(FLERR,"Group for fix_modify temp != fix group"); + return 2; + } + return 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixLangevinDrude::post_force(int /*vflag*/) +{ + // Thermalize by adding the langevin force if thermalize=true. + // Each core-Drude pair is thermalized only once: where the core is local. + + double **v = atom->v, **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal, nall = atom->nlocal + atom->nghost; + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; + double ftm2v = force->ftm2v, mvv2e = force->mvv2e; + double kb = force->boltz, dt = update->dt; + + int *drudetype = fix_drude->drudetype; + tagint *drudeid = fix_drude->drudeid; + double vdrude[3], vcore[3]; // velocities in reduced representation + double fdrude[3], fcore[3]; // forces in reduced representation + double Ccore, Cdrude, Gcore, Gdrude; + double fcoresum[3], fcoreloc[3]; + int dim = domain->dimension; + + // Compute target core temperature + if (tstyle_core == CONSTANT) + t_target_core = t_start_core; // + delta * (t_stop-t_start_core); + else { + modify->clearstep_compute(); + t_target_core = input->variable->compute_equal(tvar_core); + if (t_target_core < 0.0) + error->one(FLERR, "Fix langevin/drude variable returned " + "negative core temperature"); + modify->addstep_compute(update->ntimestep + nevery); + } + + // Compute target drude temperature + if (tstyle_drude == CONSTANT) + t_target_drude = t_start_drude; // + delta * (t_stop-t_start_core); + else { + modify->clearstep_compute(); + t_target_drude = input->variable->compute_equal(tvar_drude); + if (t_target_drude < 0.0) + error->one(FLERR, "Fix langevin/drude variable returned " + "negative drude temperature"); + modify->addstep_compute(update->ntimestep + nevery); + } + + // Clear ghost forces + // They have already been communicated if needed + for (int i = nlocal; i < nall; i++) { + for (int k = 0; k < dim; k++) + f[i][k] = 0.; + } + if (zero) for (int k=0; kremove_bias(i, v[i]); + for(int k = 0; k < dim; k++){ + fcore[k] = Ccore * random_core->gaussian() - Gcore * v[i][k]; + if (zero) fcoreloc[k] += fcore[k]; + f[i][k] += fcore[k]; + } + if (temperature) temperature->restore_bias(i, v[i]); + } else { + if (drudetype[type[i]] == DRUDE_TYPE) continue; // do with the core + + int j = atom->map(drudeid[i]); + double mi, mj, mtot, mu; // i is core, j is drude + if (rmass) { + mi = rmass[i]; + mj = rmass[j]; + } else { + mi = mass[type[i]]; + mj = mass[type[j]]; + } + mtot = mi + mj; + mu = mi * mj / mtot; + mi /= mtot; + mj /= mtot; + + Gcore = mtot / t_period_core / ftm2v; + Gdrude = mu / t_period_drude / ftm2v; + Ccore = sqrt(2.0 * Gcore * kb * t_target_core / dt / ftm2v / mvv2e); + Cdrude = sqrt(2.0 * Gdrude * kb * t_target_drude / dt / ftm2v / mvv2e); + + if (temperature) { + temperature->remove_bias(i, v[i]); + temperature->remove_bias(j, v[j]); + } + for (int k=0; kgaussian() - Gcore * vcore[k]; + fdrude[k] = Cdrude * random_drude->gaussian() - Gdrude * vdrude[k]; + + if (zero) fcoreloc[k] += fcore[k]; + + f[i][k] += mi * fcore[k] - fdrude[k]; + f[j][k] += mj * fcore[k] + fdrude[k]; + + // TODO tally energy if asked + } + if (temperature) { + temperature->restore_bias(i, v[i]); + temperature->restore_bias(j, v[j]); + } + } + } + } + + if(zero) { // Remove the drift + MPI_Allreduce(fcoreloc, fcoresum, dim, MPI_DOUBLE, MPI_SUM, world); + for (int k=0; kmap(drudeid[i]); + double mi, mj, mtot; // i is core, j is drude + if (rmass) { + mi = rmass[i]; + mj = rmass[j]; + } else { + mi = mass[type[i]]; + mj = mass[type[j]]; + } + mtot = mi + mj; + mi /= mtot; + mj /= mtot; + for (int k=0; kreverse_comm(); +} + +/* ---------------------------------------------------------------------- */ + +void FixLangevinDrude::reset_target(double t_new) +{ + t_target_core = t_start_core = t_new; +} + +/* ---------------------------------------------------------------------- + extract thermostat properties +------------------------------------------------------------------------- */ + +void *FixLangevinDrude::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"t_target_core") == 0) { + return &t_target_core; + } else if (strcmp(str,"t_target_drude") == 0) { + return &t_target_drude; + } else error->all(FLERR, "Illegal extract string in fix langevin/drude"); + return NULL; +} + +/* ---------------------------------------------------------------------- */ +int FixLangevinDrude::pack_reverse_comm(int n, int first, double *buf) +{ + int i,m,last; + double ** f = atom->f; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = f[i][0]; + buf[m++] = f[i][1]; + buf[m++] = f[i][2]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void FixLangevinDrude::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,j,m; + double ** f = atom->f; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + f[j][0] += buf[m++]; + f[j][1] += buf[m++]; + f[j][2] += buf[m++]; + } +} + diff --git a/src/USER-DRUDE/fix_langevin_drude.h b/src/USER-DRUDE/fix_langevin_drude.h new file mode 100644 index 0000000000..d65440582c --- /dev/null +++ b/src/USER-DRUDE/fix_langevin_drude.h @@ -0,0 +1,63 @@ +/* -*- 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 FIX_CLASS + +FixStyle(langevin/drude,FixLangevinDrude) + +#else + +#ifndef LMP_FIX_LANGEVIN_DRUDE_H +#define LMP_FIX_LANGEVIN_DRUDE_H + +#include "fix.h" +#include "fix_drude.h" + +namespace LAMMPS_NS { + +class FixLangevinDrude : public Fix { + public: + FixLangevinDrude(class LAMMPS *, int, char **); + virtual ~FixLangevinDrude(); + int setmask(); + void init(); + void setup(int vflag); + virtual void post_force(int vflag); + void reset_target(double); + virtual void *extract(const char *, int &); + int pack_reverse_comm(int, int, double*); + void unpack_reverse_comm(int, int*, double*); + int modify_param(int, char **); + + protected: + double t_start_core,t_period_core,t_target_core; + double t_start_drude,t_period_drude,t_target_drude; + int tstyle_core, tstyle_drude; + int tvar_core, tvar_drude; + char *tstr_core, *tstr_drude; + double energy; + int tflag; + + class RanMars *random_core, *random_drude; + int zero; + bigint ncore; + FixDrude * fix_drude; + class Compute *temperature; + char *id_temp; +}; + +} + +#endif +#endif + diff --git a/src/USER-DRUDE/pair_lj_cut_thole_long.cpp b/src/USER-DRUDE/pair_lj_cut_thole_long.cpp new file mode 100644 index 0000000000..a31a3b29f7 --- /dev/null +++ b/src/USER-DRUDE/pair_lj_cut_thole_long.cpp @@ -0,0 +1,703 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lj_cut_thole_long.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "update.h" +#include "integrate.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define EWALD_F 1.12837917 +#define EWALD_P 9.95473818e-1 +#define B0 -0.1335096380159268 +#define B1 -2.57839507e-1 +#define B2 -1.37203639e-1 +#define B3 -8.88822059e-3 +#define B4 -5.80844129e-3 +#define B5 1.14652755e-1 + +/* ---------------------------------------------------------------------- */ + +PairLJCutTholeLong::PairLJCutTholeLong(LAMMPS *lmp) : Pair(lmp) +{ + ewaldflag = pppmflag = 1; + writedata = 1; + ftable = NULL; + qdist = 0.0; + fix_drude = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairLJCutTholeLong::~PairLJCutTholeLong() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(polar); + memory->destroy(thole); + memory->destroy(ascreen); + memory->destroy(cut_lj); + memory->destroy(cut_ljsq); + memory->destroy(scale); + memory->destroy(epsilon); + memory->destroy(sigma); + memory->destroy(lj1); + memory->destroy(lj2); + memory->destroy(lj3); + memory->destroy(lj4); + memory->destroy(offset); + } + if (ftable) free_tables(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutTholeLong::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype,itable; + double qi,qj,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair,evdwl; + double r,rsq,r2inv,forcecoul,factor_coul,forcelj,factor_lj,r6inv; + double fraction,table; + double grij,expm2,prefactor,t,erfc,u; + int *ilist,*jlist,*numneigh,**firstneigh; + double factor_f,factor_e; + int di,dj; + double dqi,dqj,dcoul,asr,exp_asr; + int di_closest; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + int *drudetype = fix_drude->drudetype; + tagint *drudeid = fix_drude->drudeid; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qi = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + if (drudetype[type[i]] != NOPOL_TYPE){ + di = atom->map(drudeid[i]); + if (di < 0) error->all(FLERR, "Drude partner not found"); + di_closest = domain->closest_image(i, di); + if (drudetype[type[i]] == CORE_TYPE) + dqi = -q[di]; + else + dqi = qi; + } + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[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; + + if (rsq < cut_coulsq) { + qj = q[j]; + r = sqrt(rsq); + + if (!ncoultablebits || rsq <= tabinnersq) { + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + u = 1. - t; + erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2; + prefactor = qqrd2e * qi*qj/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = qi*qj * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = qi*qj * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + + if (drudetype[type[i]] != NOPOL_TYPE && + drudetype[type[j]] != NOPOL_TYPE){ + if (j != di_closest){ + if (drudetype[type[j]] == CORE_TYPE){ + dj = atom->map(drudeid[j]); + dqj = -q[dj]; + } else dqj = qj; + asr = ascreen[type[i]][type[j]] * r; + exp_asr = exp(-asr); + dcoul = qqrd2e * dqi * dqj / r; + factor_f = 0.5*(2. + (exp_asr * (-2. - asr * (2. + asr)))) + - factor_coul; + if (eflag) factor_e = 0.5*(2. - (exp_asr * (2. + asr))) + - factor_coul; + forcecoul += factor_f * dcoul; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + } else forcelj = 0.0; + + fpair = (forcecoul + factor_lj*forcelj) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + ecoul = prefactor*erfc; + else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qi*qj * table; + } + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + if (drudetype[type[i]] != NOPOL_TYPE && + drudetype[type[j]] != NOPOL_TYPE && j != di_closest){ + ecoul += factor_e * dcoul; + } + } else ecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairLJCutTholeLong::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); + memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); + memory->create(scale,n+1,n+1,"pair:scale"); + memory->create(ascreen,n+1,n+1,"pair:ascreen"); + memory->create(thole,n+1,n+1,"pair:thole"); + memory->create(polar,n+1,n+1,"pair:polar"); + 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 PairLJCutTholeLong::settings(int narg, char **arg) +{ + if (narg < 2 || narg > 3) error->all(FLERR,"Illegal pair_style command"); + + thole_global = force->numeric(FLERR,arg[0]); + cut_lj_global = force->numeric(FLERR,arg[1]); + if (narg == 2) cut_coul = cut_lj_global; + else cut_coul = force->numeric(FLERR,arg[2]); + + // 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]) { + thole[i][j] = thole_global; + cut_lj[i][j] = cut_lj_global; + } + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJCutTholeLong::coeff(int narg, char **arg) +{ + if (narg < 5 || narg > 7) + 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); + + double epsilon_one = force->numeric(FLERR,arg[2]); + double sigma_one = force->numeric(FLERR,arg[3]); + double polar_one = force->numeric(FLERR,arg[4]); + double thole_one = thole_global; + if (narg >=6) thole_one = force->numeric(FLERR,arg[5]); + + double cut_lj_one = cut_lj_global; + if (narg == 7) cut_lj_one = force->numeric(FLERR,arg[6]); + + 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; + polar[i][j] = polar_one; + thole[i][j] = thole_one; + ascreen[i][j] = thole[i][j] / pow(polar[i][j], 1./3.); + cut_lj[i][j] = cut_lj_one; + scale[i][j] = 1.0; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJCutTholeLong::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style lj/cut/thole/long requires atom attribute q"); + int ifix; + for (ifix = 0; ifix < modify->nfix; ifix++) + if (strcmp(modify->fix[ifix]->style,"drude") == 0) break; + if (ifix == modify->nfix) + error->all(FLERR, "Pair style lj/cut/thole/long requires fix drude"); + fix_drude = (FixDrude *) modify->fix[ifix]; + + int irequest = neighbor->request(this,instance_me); + + cut_coulsq = cut_coul * cut_coul; + + // set rRESPA cutoffs + + cut_respa = NULL; + + // insure use of KSpace long-range solver, set g_ewald + + if (force->kspace == NULL) + error->all(FLERR,"Pair style requires a KSpace style"); + g_ewald = force->kspace->g_ewald; + + // setup force tables + + if (ncoultablebits) init_tables(cut_coul,cut_respa); +} + +/* ---------------------------------------------------------------------- + neighbor callback to inform pair style of neighbor list to use + regular or rRESPA +------------------------------------------------------------------------- */ + +void PairLJCutTholeLong::init_list(int id, NeighList *ptr) +{ + if (id == 0) list = ptr; + else if (id == 1) listinner = ptr; + else if (id == 2) listmiddle = ptr; + else if (id == 3) listouter = ptr; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLJCutTholeLong::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], + sigma[i][i],sigma[j][j]); + sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); + cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); + polar[i][j] = sqrt(polar[i][i] * polar[j][j]); + thole[i][j] = 0.5 * (thole[i][i] + thole[j][j]); + ascreen[i][j] = thole[i][j] / pow(polar[i][j], 1./3.); + } + + // include TIP4P qdist in full cutoff, qdist = 0.0 if not TIP4P + + double cut = MAX(cut_lj[i][j],cut_coul+2.0*qdist); + cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; + + 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) { + double ratio = sigma[i][j] / cut_lj[i][j]; + offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); + } else offset[i][j] = 0.0; + + cut_ljsq[j][i] = cut_ljsq[i][j]; + 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]; + polar[j][i] = polar[i][j]; + thole[j][i] = thole[i][j]; + ascreen[j][i] = ascreen[i][j]; + scale[j][i] = scale[i][j]; + + // check interior rRESPA cutoff + + if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) + error->all(FLERR,"Pair cutoff < Respa interior cutoff"); + + // compute I,J contribution to long-range tail correction + // count total # of atoms of type I and J via Allreduce + + if (tail_flag) { + 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_lj[i][j]*cut_lj[i][j]*cut_lj[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; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJCutTholeLong::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(&polar[i][j],sizeof(double),1,fp); + fwrite(&thole[i][j],sizeof(double),1,fp); + fwrite(&cut_lj[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutTholeLong::read_restart(FILE *fp) +{ + read_restart_settings(fp); + + allocate(); + + 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(&polar[i][j],sizeof(double),1,fp); + fread(&thole[i][j],sizeof(double),1,fp); + ascreen[i][j] = thole[i][j] / pow(polar[i][j], 1./3.); + fread(&cut_lj[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(&polar[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&thole[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&ascreen[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJCutTholeLong::write_restart_settings(FILE *fp) +{ + fwrite(&cut_lj_global,sizeof(double),1,fp); + fwrite(&cut_coul,sizeof(double),1,fp); + fwrite(&thole_global,sizeof(double),1,fp); + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); + fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&ncoultablebits,sizeof(int),1,fp); + fwrite(&tabinner,sizeof(double),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutTholeLong::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_lj_global,sizeof(double),1,fp); + fread(&cut_coul,sizeof(double),1,fp); + fread(&thole_global,sizeof(double),1,fp); + fread(&cut_global,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + fread(&tail_flag,sizeof(int),1,fp); + fread(&ncoultablebits,sizeof(int),1,fp); + fread(&tabinner,sizeof(double),1,fp); + } + MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); + MPI_Bcast(&thole_global,1,MPI_DOUBLE,0,world); + 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); + MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); + MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); +} + + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairLJCutTholeLong::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairLJCutTholeLong::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut_lj[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJCutTholeLong::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, + double factor_lj, double &fforce) +{ + double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor,u; + double fraction,table,forcecoul,forcelj,phicoul,philj; + int itable; + double factor_f,factor_e; + double dqi,dqj,dcoul,asr,exp_asr; + int di, dj, di_closest; + + int *drudetype = fix_drude->drudetype; + tagint *drudeid = fix_drude->drudeid; + int *type = atom->type; + + r2inv = 1.0/rsq; + if (rsq < cut_coulsq) { + r = sqrt(rsq); + if (!ncoultablebits || rsq <= tabinnersq) { + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + u = 1. - t; + erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2; + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup_single; + rsq_lookup_single.f = rsq; + itable = rsq_lookup_single.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = atom->q[i]*atom->q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = atom->q[i]*atom->q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + if (drudetype[type[i]] != NOPOL_TYPE && drudetype[type[j]] != NOPOL_TYPE) { + di = atom->map(drudeid[i]); + di_closest = domain->closest_image(i, di); + if (di_closest != dj){ + if (drudetype[i] == CORE_TYPE) dqi = -atom->q[di]; + else if (drudetype[i] == DRUDE_TYPE) dqi = atom->q[i]; + if (drudetype[j] == CORE_TYPE) { + dj = atom->map(drudeid[j]); + dqj = -atom->q[dj]; + } else if (drudetype[j] == DRUDE_TYPE) dqj = atom->q[j]; + asr = ascreen[itype][jtype] * r; + exp_asr = exp(-asr); + dcoul = force->qqrd2e * dqi * dqj / r; + factor_f = 0.5*(2. + (exp_asr * (-2. - asr * (2. + asr)))) + - factor_coul; + forcecoul += factor_f * dcoul; + factor_e = 0.5*(2. - (exp_asr * (2. + asr))) - factor_coul; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + } else forcelj = 0.0; + + fforce = (forcecoul + factor_lj*forcelj) * r2inv; + + double eng = 0.0; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*erfc; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = atom->q[i]*atom->q[j] * table; + } + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + if (drudetype[type[i]] != NOPOL_TYPE && drudetype[type[j]] != NOPOL_TYPE && + di_closest != dj) + phicoul += factor_e * dcoul; + eng += phicoul; + } + + if (rsq < cut_ljsq[itype][jtype]) { + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + eng += factor_lj*philj; + } + + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJCutTholeLong::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + dim = 6; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str,"scale") == 0) return (void *) scale; + if (strcmp(str,"polar") == 0) return (void *) polar; + if (strcmp(str,"thole") == 0) return (void *) thole; + if (strcmp(str,"ascreen") == 0) return (void *) ascreen; + return NULL; +} diff --git a/src/USER-DRUDE/pair_lj_cut_thole_long.h b/src/USER-DRUDE/pair_lj_cut_thole_long.h new file mode 100644 index 0000000000..894042f6ce --- /dev/null +++ b/src/USER-DRUDE/pair_lj_cut_thole_long.h @@ -0,0 +1,97 @@ +/* -*- 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 PAIR_CLASS + +PairStyle(lj/cut/thole/long,PairLJCutTholeLong) + +#else + +#ifndef LMP_PAIR_LJ_CUT_THOLE_LONG_H +#define LMP_PAIR_LJ_CUT_THOLE_LONG_H + +#include "pair.h" +#include "fix_drude.h" + +namespace LAMMPS_NS { + +class PairLJCutTholeLong : public Pair { + + public: + PairLJCutTholeLong(class LAMMPS *); + virtual ~PairLJCutTholeLong(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + virtual void init_style(); + void init_list(int, class NeighList *); + virtual double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + virtual void write_restart_settings(FILE *); + virtual void read_restart_settings(FILE *); + void write_data(FILE *); + void write_data_all(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); + + virtual void *extract(const char *, int &); + + protected: + double cut_lj_global; + double **cut_lj,**cut_ljsq; + double cut_coul,cut_coulsq; + double **epsilon,**sigma; + double **lj1,**lj2,**lj3,**lj4,**offset; + double *cut_respa; + double qdist; // TIP4P distance from O site to negative charge + double g_ewald; + double thole_global; + double cut_global; + double **cut,**scale; + double **polar,**thole,**ascreen; + FixDrude *fix_drude; + + virtual 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: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style lj/cut/coul/long requires atom attribute q + +The atom style defined does not have this attribute. + +E: Pair style requires a KSpace style + +No kspace style is defined. + +E: Pair cutoff < Respa interior cutoff + +One or more pairwise cutoffs are too short to use with the specified +rRESPA cutoffs. + +*/ diff --git a/src/USER-DRUDE/pair_thole.cpp b/src/USER-DRUDE/pair_thole.cpp new file mode 100644 index 0000000000..90e18f8742 --- /dev/null +++ b/src/USER-DRUDE/pair_thole.cpp @@ -0,0 +1,424 @@ +/* ---------------------------------------------------------------------- + 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 "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_thole.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" +#include "fix.h" +#include "fix_store.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairThole::PairThole(LAMMPS *lmp) : Pair(lmp) { + fix_drude = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairThole::~PairThole() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(polar); + memory->destroy(thole); + memory->destroy(ascreen); + memory->destroy(cut); + memory->destroy(scale); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairThole::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qi,qj,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; + double r,rsq,r2inv,rinv,forcecoul,factor_coul; + int *ilist,*jlist,*numneigh,**firstneigh; + double factor_f,factor_e; + int di,dj; + double dcoul,asr,exp_asr; + + ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + int *drudetype = fix_drude->drudetype; + tagint *drudeid = fix_drude->drudeid; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + + // only on core-drude pair + if (drudetype[type[i]] == NOPOL_TYPE) + continue; + + di = domain->closest_image(i, atom->map(drudeid[i])); + // get dq of the core via the drude charge + if (drudetype[type[i]] == DRUDE_TYPE) + qi = q[i]; + else + qi = -q[di]; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + // only on core-drude pair, but not into the same pair + if (drudetype[type[j]] == NOPOL_TYPE || j == di) + continue; + + // get dq of the core via the drude charge + if (drudetype[type[j]] == DRUDE_TYPE) + qj = q[j]; + else { + dj = domain->closest_image(j, atom->map(drudeid[j])); + qj = -q[dj]; + } + + 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; + rinv = sqrt(r2inv); + + r = sqrt(rsq); + asr = ascreen[itype][jtype] * r; + exp_asr = exp(-asr); + dcoul = qqrd2e * qi * qj *scale[itype][jtype] * rinv; + factor_f = 0.5*(2. + (exp_asr * (-2. - asr * (2. + asr)))) - factor_coul; + if(eflag) factor_e = 0.5*(2. - (exp_asr * (2. + asr))) - factor_coul; + fpair = factor_f * dcoul * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) + ecoul = factor_e * dcoul; + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + 0.0,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairThole::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(scale,n+1,n+1,"pair:scale"); + memory->create(ascreen,n+1,n+1,"pair:ascreen"); + memory->create(thole,n+1,n+1,"pair:thole"); + memory->create(polar,n+1,n+1,"pair:polar"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairThole::settings(int narg, char **arg) +{ + if (narg != 2) error->all(FLERR,"Illegal pair_style command"); + + thole_global = force->numeric(FLERR,arg[0]); + cut_global = force->numeric(FLERR,arg[1]); + + // 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]) { + thole[i][j] = thole_global; + cut[i][j] = cut_global; + } + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairThole::coeff(int narg, char **arg) +{ + if (narg < 3 || narg > 5) + 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); + + double polar_one = force->numeric(FLERR,arg[2]); + double thole_one = thole_global; + double cut_one = cut_global; + if (narg >=4) thole_one = force->numeric(FLERR,arg[3]); + if (narg == 5) cut_one = force->numeric(FLERR,arg[4]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + polar[i][j] = polar_one; + thole[i][j] = thole_one; + ascreen[i][j] = thole[i][j] / pow(polar[i][j], 1./3.); + cut[i][j] = cut_one; + scale[i][j] = 1.0; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairThole::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style thole requires atom attribute q"); + int ifix; + for (ifix = 0; ifix < modify->nfix; ifix++) + if (strcmp(modify->fix[ifix]->style,"drude") == 0) break; + if (ifix == modify->nfix) error->all(FLERR, "Pair thole requires fix drude"); + fix_drude = (FixDrude *) modify->fix[ifix]; + + neighbor->request(this,instance_me); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairThole::init_one(int i, int j) +{ + if (setflag[i][j] == 0) + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + + polar[j][i] = polar[i][j]; + thole[j][i] = thole[i][j]; + ascreen[j][i] = ascreen[i][j]; + scale[j][i] = scale[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairThole::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(&polar[i][j],sizeof(double),1,fp); + fwrite(&thole[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairThole::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + 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(&polar[i][j],sizeof(double),1,fp); + fread(&thole[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + ascreen[i][j] = thole[i][j] / pow(polar[i][j], 1./3.); + } + MPI_Bcast(&polar[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&thole[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&ascreen[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairThole::write_restart_settings(FILE *fp) +{ + fwrite(&thole_global,sizeof(double),1,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 PairThole::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&thole_global,sizeof(double),1,fp); + fread(&cut_global,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&thole_global,1,MPI_DOUBLE,0,world); + 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 PairThole::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,rinv,r,forcecoul,phicoul; + double qi,qj,factor_f,factor_e,dcoul,asr,exp_asr; + int di, dj; + + int *drudetype = fix_drude->drudetype; + tagint *drudeid = fix_drude->drudeid; + int *type = atom->type; + + // only on core-drude pair, but not on the same pair + if (drudetype[type[i]] == NOPOL_TYPE || drudetype[type[j]] == NOPOL_TYPE || + j == i) + return 0.0; + + // get dq of the core via the drude charge + if (drudetype[type[i]] == DRUDE_TYPE) + qi = atom->q[i]; + else { + di = domain->closest_image(i, atom->map(drudeid[i])); + qi = -atom->q[di]; + } + if (drudetype[type[j]] == DRUDE_TYPE) + qj = atom->q[j]; + else { + dj = domain->closest_image(j, atom->map(drudeid[j])); + qj = -atom->q[dj]; + } + + r2inv = 1.0/rsq; + if (rsq < cutsq[itype][jtype]) { + rinv = sqrt(r2inv); + r = sqrt(rsq); + asr = ascreen[itype][jtype] * r; + exp_asr = exp(-asr); + dcoul = force->qqrd2e * qi * qj * rinv; + factor_f = 0.5*(2. + (exp_asr * (-2. - asr * (2. + asr)))) - factor_coul; + forcecoul += factor_f * dcoul; + factor_e = 0.5*(2. - (exp_asr * (2. + asr))) - factor_coul; + } else forcecoul= 0.0; + fforce = factor_f*forcecoul * r2inv; + + if (rsq < cutsq[itype][jtype]) { + phicoul = factor_e * dcoul; + } else phicoul = 0.0; + + return phicoul; +} + +/* ---------------------------------------------------------------------- */ + +void *PairThole::extract(const char *str, int &dim) +{ + dim = 4; + if (strcmp(str,"scale") == 0) return (void *) scale; + if (strcmp(str,"polar") == 0) return (void *) polar; + if (strcmp(str,"thole") == 0) return (void *) thole; + if (strcmp(str,"ascreen") == 0) return (void *) ascreen; + return NULL; +} diff --git a/src/USER-DRUDE/pair_thole.h b/src/USER-DRUDE/pair_thole.h new file mode 100644 index 0000000000..2c462b431a --- /dev/null +++ b/src/USER-DRUDE/pair_thole.h @@ -0,0 +1,75 @@ +/* -*- 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 PAIR_CLASS + +PairStyle(thole,PairThole) + +#else + +#ifndef LMP_PAIR_THOLE_H +#define LMP_PAIR_THOLE_H + +#include "pair.h" +#include "fix_drude.h" + +namespace LAMMPS_NS { + +class PairThole : public Pair { + public: + PairThole(class LAMMPS *); + virtual ~PairThole(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + void init_style(); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + virtual void write_restart_settings(FILE *); + virtual void read_restart_settings(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); + void *extract(const char *, int &); + + protected: + double thole_global; + double cut_global; + double **cut,**scale; + double **polar,**thole,**ascreen; + FixDrude * fix_drude; + + virtual 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: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style thole requires atom attribute q + +The atom style defined does not have these attributes. + +*/