From 8889d8dfb96dfe979bbc620aee1196676aceedb1 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 22 Oct 2010 21:26:02 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5103 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_ti.cpp | 180 +++++++++++++++++++++++++++++++++++++++++++++ src/compute_ti.h | 46 ++++++++++++ src/pair_born.cpp | 11 +++ src/pair_born.h | 1 + src/pair_buck.cpp | 11 +++ src/pair_buck.h | 1 + src/pair_gauss.cpp | 32 ++++++-- src/pair_gauss.h | 42 ++++++----- 8 files changed, 299 insertions(+), 25 deletions(-) create mode 100644 src/compute_ti.cpp create mode 100644 src/compute_ti.h diff --git a/src/compute_ti.cpp b/src/compute_ti.cpp new file mode 100644 index 0000000000..fb9ed99359 --- /dev/null +++ b/src/compute_ti.cpp @@ -0,0 +1,180 @@ +/* ---------------------------------------------------------------------- + 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: Sai Jayaraman (University of Notre Dame) +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "string.h" +#include "compute_ti.h" +#include "update.h" +#include "modify.h" +#include "domain.h" +#include "force.h" +#include "pair.h" +#include "kspace.h" +#include "input.h" +#include "variable.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{PAIR,TAIL,KSPACE}; + +#define INVOKED_SCALAR 1 + +/* ---------------------------------------------------------------------- */ + +ComputeTI::ComputeTI(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 4) error->all("Illegal compute ti command"); + + peflag = 1; + scalar_flag = 1; + extscalar = 1; + timeflag = 1; + + // terms come in triplets + + nterms = (narg-3) / 3; + if (narg != 3*nterms + 3) error->all("Illegal compute ti command"); + + which = new int[nterms]; + ivar1 = new int[nterms]; + ivar2 = new int[nterms]; + var1 = new char*[nterms]; + var2 = new char*[nterms]; + pptr = new Pair*[nterms]; + pstyle = new char*[nterms]; + + for (int m = 0; m < nterms; m++) pstyle[m] = NULL; + + // parse keywords + + nterms = 0; + + int iarg = 3; + while (iarg < narg) { + if (iarg+3 > narg) error->all("Illegal fix adapt command"); + if (strcmp(arg[iarg],"kspace") == 0) which[nterms] = KSPACE; + else if (strcmp(arg[iarg],"tail") == 0) which[nterms] = TAIL; + else { + which[nterms] = PAIR; + int n = strlen(&arg[iarg+1][2]) + 1; + pstyle[nterms] = new char[n]; + strcpy(pstyle[nterms],&arg[iarg+1][2]); + } + + if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { + int n = strlen(&arg[iarg+2][2]) + 1; + var1[nterms] = new char[n]; + strcpy(var1[nterms],&arg[iarg+2][2]); + } else error->all("Illegal compute ti command"); + if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { + int n = strlen(&arg[iarg+3][2]) + 1; + var2[nterms] = new char[n]; + strcpy(var2[nterms],&arg[iarg+3][2]); + } else error->all("Illegal compute ti command"); + + nterms++; + iarg += 3; + } +} + +/* --------------------------------------------------------------------- */ + +ComputeTI::~ComputeTI() +{ + for (int m = 0; m < nterms; m++) { + delete [] var1[m]; + delete [] var2[m]; + delete [] pstyle[m]; + } + delete [] which; + delete [] ivar1; + delete [] ivar2; + delete [] var1; + delete [] var2; + delete [] pptr; + delete [] pstyle; +} + +/* --------------------------------------------------------------------- */ + +void ComputeTI::init() +{ + // setup and error checks + + for (int m = 0; m < nterms; m++) { + ivar1[m] = input->variable->find(var1[m]); + ivar2[m] = input->variable->find(var2[m]); + if (ivar1[m] < 0 || ivar2 < 0) + error->all("Variable name for compute ti does not exist"); + if (!input->variable->equalstyle(ivar1[m]) || + !input->variable->equalstyle(ivar2[m])) + error->all("Variable for compute ti is invalid style"); + + if (which[m] == PAIR) { + pptr[m] = force->pair_match(pstyle[m],1); + if (pptr[m] == NULL) error->all("Compute ti pair style does not exist"); + + } else if (which[m] == TAIL) { + if (force->pair == NULL || force->pair->tail_flag == 0) + error->all("Compute ti tail when pair style does not " + "compute tail corrections"); + + } else if (which[m] == KSPACE) { + if (force->kspace == NULL) + error->all("Compute ti is incompatible with KSpace style"); + } + } +} + +/* --------------------------------------------------------------------- */ + +double ComputeTI::compute_scalar() +{ + double eng,engall,value1,value2; + + invoked_scalar |= INVOKED_SCALAR; + if (update->eflag_global != invoked_scalar) + error->all("Energy was not tallied on needed timestep"); + + double dUdl = 0.0; + + for (int m = 0; m < nterms; m++) { + value1 = input->variable->compute_equal(ivar1[m]); + value2 = input->variable->compute_equal(ivar2[m]); + if (value1 == 0.0) continue; + + if (which[m] == PAIR) { + eng = pptr[m]->eng_vdwl + pptr[m]->eng_coul; + MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); + dUdl += engall/value1 * value2; + + } else if (which[m] == TAIL) { + double vol = domain->xprd*domain->yprd*domain->zprd; + eng = force->pair->etail / vol; + dUdl += eng/value1 * value2; + + } else if (which[m] == KSPACE) { + eng = force->kspace->energy; + dUdl += eng/value1 * value2; + } + } + + scalar = dUdl; + return scalar; +} diff --git a/src/compute_ti.h b/src/compute_ti.h new file mode 100644 index 0000000000..02dfacff17 --- /dev/null +++ b/src/compute_ti.h @@ -0,0 +1,46 @@ +/* ---------------------------------------------------------------------- + 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(ti,ComputeTI) + +#else + +#ifndef COMPUTE_TI_H +#define COMPUTE_TI_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeTI : public Compute { + public: + ComputeTI(class LAMMPS *, int, char **); + ~ComputeTI(); + void init(); + double compute_scalar(); + + private: + int nterms; + int *which; + int *ivar1,*ivar2; + char **var1,**var2; + class Pair **pptr; + char **pstyle; +}; + +} + +#endif +#endif diff --git a/src/pair_born.cpp b/src/pair_born.cpp index f015552252..13194d43b0 100644 --- a/src/pair_born.cpp +++ b/src/pair_born.cpp @@ -403,3 +403,14 @@ double PairBorn::single(int i, int j, int itype, int jtype, d[itype][jtype]*r2inv*r6inv - offset[itype][jtype]; return factor_lj*phiborn; } + +/* ---------------------------------------------------------------------- */ + +void *PairBorn::extract(char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"a") == 0) return (void *) a; + if (strcmp(str,"c") == 0) return (void *) c; + if (strcmp(str,"d") == 0) return (void *) d; + return NULL; +} diff --git a/src/pair_born.h b/src/pair_born.h index 8b62e47f12..20fe5db8cb 100644 --- a/src/pair_born.h +++ b/src/pair_born.h @@ -37,6 +37,7 @@ class PairBorn : public Pair { void write_restart_settings(FILE *); void read_restart_settings(FILE *); double single(int, int, int, int, double, double, double, double &); + void *extract(char *, int &); private: double cut_global; diff --git a/src/pair_buck.cpp b/src/pair_buck.cpp index 875a749d47..f8c48ca1c5 100644 --- a/src/pair_buck.cpp +++ b/src/pair_buck.cpp @@ -14,6 +14,7 @@ #include "math.h" #include "stdio.h" #include "stdlib.h" +#include "string.h" #include "pair_buck.h" #include "atom.h" #include "comm.h" @@ -342,3 +343,13 @@ double PairBuck::single(int i, int j, int itype, int jtype, offset[itype][jtype]; return factor_lj*phibuck; } + +/* ---------------------------------------------------------------------- */ + +void *PairBuck::extract(char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"a") == 0) return (void *) a; + if (strcmp(str,"c") == 0) return (void *) c; + return NULL; +} diff --git a/src/pair_buck.h b/src/pair_buck.h index 21ba0563f9..70d571c8b8 100644 --- a/src/pair_buck.h +++ b/src/pair_buck.h @@ -37,6 +37,7 @@ class PairBuck : public Pair { void write_restart_settings(FILE *); void read_restart_settings(FILE *); double single(int, int, int, int, double, double, double, double &); + void *extract(char *, int &); private: double cut_global; diff --git a/src/pair_gauss.cpp b/src/pair_gauss.cpp index c9607d815c..f59c3ad270 100644 --- a/src/pair_gauss.cpp +++ b/src/pair_gauss.cpp @@ -37,12 +37,18 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairGauss::PairGauss(LAMMPS *lmp) :Pair(lmp) {} +PairGauss::PairGauss(LAMMPS *lmp) :Pair(lmp) +{ + nextra = 1; + pvector = new double[1]; +} /* ---------------------------------------------------------------------- */ PairGauss::~PairGauss() { + delete [] pvector; + if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); @@ -66,7 +72,8 @@ void PairGauss::compute(int eflag, int vflag) evdwl = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; - + int occ = 0; + double **x = atom->x; double **f = atom->f; int *type = atom->type; @@ -98,7 +105,12 @@ void PairGauss::compute(int eflag, int vflag) delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; - + + // define a Gaussian well to be occupied if + // the site it interacts with is within the force maximum + + if (eflag_global && rsq < 0.5/b[itype][jtype]) occ++; + if (rsq < cutsq[itype][jtype]) { r2inv = 1.0/rsq; r = sqrt(rsq); @@ -120,11 +132,12 @@ void PairGauss::compute(int eflag, int vflag) offset[itype][jtype]); if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); + evdwl,0.0,fpair,delx,dely,delz); } } } - + + if (eflag_global) pvector[0] = occ; if (vflag_fdotr) virial_compute(); } @@ -314,3 +327,12 @@ double PairGauss::single(int i, int j, int itype, int jtype, double rsq, fforce = forcelj*r2inv; return philj; } + +/* ---------------------------------------------------------------------- */ + +void *PairGauss::extract(char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"a") == 0) return (void *) a; + return NULL; +} diff --git a/src/pair_gauss.h b/src/pair_gauss.h index 9598457dca..c8e94a111b 100644 --- a/src/pair_gauss.h +++ b/src/pair_gauss.h @@ -24,28 +24,30 @@ PairStyle(gauss,PairGauss) namespace LAMMPS_NS { - class PairGauss : public Pair { - public: - PairGauss(class LAMMPS *); - ~PairGauss(); - void compute(int,int); - void settings(int, char **); - void coeff(int, char **); - 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 &); +class PairGauss : public Pair { + public: + PairGauss(class LAMMPS *); + ~PairGauss(); + void compute(int,int); + void settings(int, char **); + void coeff(int, char **); + 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 &); - private: - double cut_global; - double **cut; - double **a,**b; - double **offset; + private: + double cut_global; + double **cut; + double **a,**b; + double **offset; + + void allocate(); +}; - void allocate(); - }; } #endif