diff --git a/src/Makefile b/src/Makefile index 019bb4bed9..8241135cc2 100755 --- a/src/Makefile +++ b/src/Makefile @@ -18,7 +18,7 @@ PACKAGE = asphere body class2 colloid dipole fld gpu granular kim \ reax replica rigid shock srd voronoi xtc PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \ - user-cuda user-eff user-lb user-misc user-molfile \ + user-cuda user-eff user-fep user-lb user-misc user-molfile \ user-omp user-phonon user-qmmm user-reaxc user-sph PACKLIB = gpu kim meam poems reax voronoi \ diff --git a/src/USER-FEP/Install.sh b/src/USER-FEP/Install.sh new file mode 100644 index 0000000000..b89adacd36 --- /dev/null +++ b/src/USER-FEP/Install.sh @@ -0,0 +1,48 @@ +# Install/unInstall package files in LAMMPS +# mode = 0/1/2 for uninstall/install/update + +# this is default Install.sh for all packages +# if package has an auxiliary library or a file with a dependency, +# then package dir has its own customized Install.sh + +mode=$1 + +# arg1 = file, arg2 = file it depends on + +action () { + if (test $mode = 0) then + rm -f ../$1 + elif (! cmp -s $1 ../$1) then + if (test -z "$2" || test -e ../$2) then + cp $1 .. + if (test $mode = 2) then + echo " updating src/$1" + fi + fi + elif (test -n "$2") then + if (test ! -e ../$2) then + rm -f ../$1 + fi + fi +} + +# all package files with dependencies + +action compute_fep.cpp +action compute_fep.h +action fix_adapt_fep.cpp +action fix_adapt_fep.h +action pair_coul_cut_soft.cpp +action pair_coul_cut_soft.h +action pair_coul_long_soft.cpp pppm.cpp +action pair_coul_long_soft.h pppm.cpp +action pair_lj_cut_coul_cut_soft.cpp +action pair_lj_cut_coul_cut_soft.h +action pair_lj_cut_coul_long_soft.cpp pppm.cpp +action pair_lj_cut_coul_long_soft.h pppm.cpp +action pair_lj_cut_soft.cpp +action pair_lj_cut_soft.h +action pair_lj_cut_tip4p_long_soft.cpp pppm_tip4p.cpp +action pair_lj_cut_tip4p_long_soft.h pppm_tip4p.cpp +action pair_tip4p_long_soft.cpp pppm_tip4p.cpp +action pair_tip4p_long_soft.h pppm_tip4p.cpp diff --git a/src/USER-FEP/README b/src/USER-FEP/README new file mode 100644 index 0000000000..967f8e6efe --- /dev/null +++ b/src/USER-FEP/README @@ -0,0 +1,13 @@ +This package provides methods for performing free energy perturbation +simulations with soft-core pair potentials in LAMMPS. + +See the doc page for the fix adapt/fep command to get started. + +There are example scripts for using this package in +examples/USER/fep. + +There are auxiliary tools for using this package in tools/fep. + +The person who created this package is Agilio Padua at Université +Blaise Pascal Clermont-Ferrand (agilio.padua at univ-bpclermont.fr) +Contact him directly if you have questions. diff --git a/src/USER-FEP/compute_fep.cpp b/src/USER-FEP/compute_fep.cpp new file mode 100644 index 0000000000..847502cabc --- /dev/null +++ b/src/USER-FEP/compute_fep.cpp @@ -0,0 +1,690 @@ +/* ---------------------------------------------------------------------- + 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: Agilio Padua (Univ Blaise Pascal & CNRS) +------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "string.h" +#include "math.h" +#include "mpi.h" +#include "comm.h" +#include "update.h" +#include "atom.h" +#include "domain.h" +#include "force.h" +#include "pair.h" +#include "pair_hybrid.h" +#include "kspace.h" +#include "input.h" +#include "variable.h" +#include "timer.h" +#include "memory.h" +#include "error.h" +#include "compute_fep.h" + +using namespace LAMMPS_NS; + +enum{PAIR,ATOM}; +enum{CHARGE}; + +#undef FEP_DEBUG +#undef FEP_MAXDEBUG + +/* ---------------------------------------------------------------------- */ + +ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 5) error->all(FLERR,"Illegal number of arguments in compute fep"); + + scalar_flag = 0; + vector_flag = 1; + size_vector = 3; + extvector = 0; + + vector = new double[3]; + + fepinitflag = 0; // avoid init to run entirely when called by write_data + + temp_fep = force->numeric(FLERR,arg[3]); + + // count # of perturbations + + npert = 0; + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"pair") == 0) { + if (iarg+6 > narg) error->all(FLERR,"Illegal pair attribute in compute fep"); + npert++; + iarg += 6; + } else if (strcmp(arg[iarg],"atom") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal atom attribute in compute fep"); + npert++; + iarg += 4; + } else break; + } + + if (npert == 0) error->all(FLERR,"Illegal syntax in compute fep"); + perturb = new Perturb[npert]; + + // parse keywords + + npert = 0; + chgflag = 0; + + iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"pair") == 0) { + perturb[npert].which = PAIR; + int n = strlen(arg[iarg+1]) + 1; + perturb[npert].pstyle = new char[n]; + strcpy(perturb[npert].pstyle,arg[iarg+1]); + n = strlen(arg[iarg+2]) + 1; + perturb[npert].pparam = new char[n]; + strcpy(perturb[npert].pparam,arg[iarg+2]); + force->bounds(arg[iarg+3],atom->ntypes, + perturb[npert].ilo,perturb[npert].ihi); + force->bounds(arg[iarg+4],atom->ntypes, + perturb[npert].jlo,perturb[npert].jhi); + if (strstr(arg[iarg+5],"v_") == arg[iarg+5]) { + n = strlen(&arg[iarg+5][2]) + 1; + perturb[npert].var = new char[n]; + strcpy(perturb[npert].var,&arg[iarg+5][2]); + } else error->all(FLERR,"Illegal variable in compute fep"); + npert++; + iarg += 6; + } else if (strcmp(arg[iarg],"atom") == 0) { + perturb[npert].which = ATOM; + if (strcmp(arg[iarg+1],"charge") == 0) { + perturb[npert].aparam = CHARGE; + chgflag = 1; + } else error->all(FLERR,"Illegal atom argument in compute fep"); + force->bounds(arg[iarg+2],atom->ntypes, + perturb[npert].ilo,perturb[npert].ihi); + if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { + int n = strlen(&arg[iarg+3][2]) + 1; + perturb[npert].var = new char[n]; + strcpy(perturb[npert].var,&arg[iarg+3][2]); + } else error->all(FLERR,"Illegal variable in compute fep"); + npert++; + iarg += 4; + } else break; + } + + // optional keywords + + tailflag = 0; + volumeflag = 0; + + while (iarg < narg) { + if (strcmp(arg[iarg],"tail") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal optional keyword " + "in compute fep"); + if (strcmp(arg[iarg+1],"no") == 0) tailflag = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) tailflag = 1; + else error->all(FLERR,"Illegal optional keyword in compute fep"); + iarg += 2; + } else if (strcmp(arg[iarg],"volume") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal optional keyword " + "in compute fep"); + if (strcmp(arg[iarg+1],"no") == 0) volumeflag = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) volumeflag = 1; + else error->all(FLERR,"Illegal optional keyword in compute fep"); + iarg += 2; + } else + error->all(FLERR,"Illegal optional keyword in compute fep"); + } + + // allocate pair style arrays + + int ntype = atom->ntypes; + for (int m = 0; m < npert; m++) { + if (perturb[m].which == PAIR) + memory->create(perturb[m].array_orig,ntype+1,ntype+1,"fep:array_orig"); + } + + // allocate space for charge, force, energy, virial arrays + + allocate_storage(); + +} + +/* ---------------------------------------------------------------------- */ + +ComputeFEP::~ComputeFEP() +{ + delete [] vector; + + for (int m = 0; m < npert; m++) { + delete [] perturb[m].var; + if (perturb[m].which == PAIR) { + delete [] perturb[m].pstyle; + delete [] perturb[m].pparam; + memory->destroy(perturb[m].array_orig); + } + } + delete [] perturb; + + deallocate_storage(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFEP::init() +{ + int i,j; + + if (!fepinitflag) // avoid init to run entirely when called by write_data + fepinitflag = 1; + else return; + + // setup and error checks + + pairflag = 0; + + for (int m = 0; m < npert; m++) { + Perturb *pert = &perturb[m]; + + pert->ivar = input->variable->find(pert->var); + if (pert->ivar < 0) + error->all(FLERR,"Variable name for compute fep does not exist"); + if (!input->variable->equalstyle(pert->ivar)) + error->all(FLERR,"Variable for compute fep is of invalid style"); + + if (force->pair == NULL) + error->all(FLERR,"compute fep pair requires pair interactions"); + + if (pert->which == PAIR) { + pairflag = 1; + + Pair *pair = force->pair_match(pert->pstyle,1); + if (pair == NULL) error->all(FLERR,"compute fep pair style " + "does not exist"); + void *ptr = pair->extract(pert->pparam,pert->pdim); + if (ptr == NULL) + error->all(FLERR,"compute fep pair style param not supported"); + + pert->array = (double **) ptr; + + // if pair hybrid, test that ilo,ihi,jlo,jhi are valid for sub-style + + if ((strcmp(force->pair_style,"hybrid") == 0 || + strcmp(force->pair_style,"hybrid/overlay") == 0)) { + PairHybrid *pair = (PairHybrid *) force->pair; + for (i = pert->ilo; i <= pert->ihi; i++) + for (j = MAX(pert->jlo,i); j <= pert->jhi; j++) + if (!pair->check_ijtype(i,j,pert->pstyle)) + error->all(FLERR,"compute fep type pair range is not valid for " + "pair hybrid sub-style"); + } + + } else if (pert->which == ATOM) { + if (pert->aparam == CHARGE) { + if (!atom->q_flag) + error->all(FLERR,"compute fep requires atom attribute charge"); + } + } + } + + if (tailflag) { + if (force->pair->tail_flag == 0) + error->all(FLERR,"Compute fep tail when pair style does not " + "compute tail corrections"); + } + + if (comm->me == 0) { + if (screen) { + fprintf(screen, "FEP settings ...\n"); + fprintf(screen, " temperature = %f\n", temp_fep); + fprintf(screen, " tail %s\n", (tailflag ? "yes":"no")); + for (int m = 0; m < npert; m++) { + Perturb *pert = &perturb[m]; + if (pert->which == PAIR) + fprintf(screen, " %s %s %d-%d %d-%d\n", pert->pstyle, pert->pparam, + pert->ilo, pert->ihi, pert->jlo, pert->jhi); + else if (pert->which == ATOM) + fprintf(screen, " %d-%d charge\n", pert->ilo, pert->ihi); + } + } + if (logfile) { + fprintf(logfile, "FEP settings ...\n"); + fprintf(logfile, " temperature = %f\n", temp_fep); + fprintf(logfile, " tail %s\n", (tailflag ? "yes":"no")); + for (int m = 0; m < npert; m++) { + Perturb *pert = &perturb[m]; + if (pert->which == PAIR) + fprintf(logfile, " %s %s %d-%d %d-%d\n", pert->pstyle, pert->pparam, + pert->ilo, pert->ihi, pert->jlo, pert->jhi); + else if (pert->which == ATOM) + fprintf(logfile, " %d-%d charge\n", pert->ilo, pert->ihi); + } + } + } + +} + +/* ---------------------------------------------------------------------- */ + + +void ComputeFEP::compute_vector() +{ + double pe0,pe1; + + invoked_vector = update->ntimestep; + + if (atom->nmax > nmax) { // reallocate working arrays if necessary + deallocate_storage(); + allocate_storage(); + } + + backup_qfev(); // backup charge, force, energy, virial array values + backup_params(); // backup pair parameters + + timer->stamp(); + if (force->pair && force->pair->compute_flag) { + force->pair->compute(1,0); + timer->stamp(TIME_PAIR); + } + if (force->kspace && force->kspace->compute_flag) { + force->kspace->compute(1,0); + timer->stamp(TIME_KSPACE); + } + pe0 = compute_epair(); + + perturb_params(); + + timer->stamp(); + if (force->pair && force->pair->compute_flag) { + force->pair->compute(1,0); + timer->stamp(TIME_PAIR); + } + if (force->kspace && force->kspace->compute_flag) { + force->kspace->compute(1,0); + timer->stamp(TIME_KSPACE); + } + pe1 = compute_epair(); + + restore_qfev(); // restore charge, force, energy, virial array values + restore_params(); // restore pair parameters + + vector[0] = pe1-pe0; + vector[1] = exp(-(pe1-pe0)/(force->boltz*temp_fep)); + vector[2] = domain->xprd * domain->yprd * domain->zprd; + if (volumeflag) + vector[1] *= vector[2]; + +#ifdef FEP_DEBUG + if (comm->me == 0 && screen) + fprintf(screen, "FEP U0 = %f U1 = %f DU = %f exp(-DU/kT) = %f\n", + pe0,pe1,vector[0],vector[1]); +#endif +} + + +/* ---------------------------------------------------------------------- + obtain pair energy from lammps accumulators +------------------------------------------------------------------------- */ + +double ComputeFEP::compute_epair() +{ + double eng, eng_pair; + + eng = 0.0; + if (force->pair) + eng = force->pair->eng_vdwl + force->pair->eng_coul; + MPI_Allreduce(&eng,&eng_pair,1,MPI_DOUBLE,MPI_SUM,world); + + if (tailflag) { + double volume = domain->xprd * domain->yprd * domain->zprd; + eng_pair += force->pair->etail / volume; + } + + if (force->kspace) eng_pair += force->kspace->energy; + + return eng_pair; +} + + +/* ---------------------------------------------------------------------- + apply perturbation to pair, atom parameters based on variable evaluation +------------------------------------------------------------------------- */ + +void ComputeFEP::perturb_params() +{ + int i,j; + + for (int m = 0; m < npert; m++) { + Perturb *pert = &perturb[m]; + + double delta = input->variable->compute_equal(pert->ivar); + + if (pert->which == PAIR) { // modify pair parameters + for (i = pert->ilo; i <= pert->ihi; i++) + for (j = MAX(pert->jlo,i); j <= pert->jhi; j++) + pert->array[i][j] = pert->array_orig[i][j] + delta; + +#ifdef FEP_MAXDEBUG + if (comm->me == 0 && screen) { + fprintf(screen, "###FEP change %s %s, delta = %f\n", + pert->pstyle, pert->pparam, delta); + fprintf(screen, "###FEP I J old_param new_param\n"); + for (i = pert->ilo; i <= pert->ihi; i++) + for (j = MAX(pert->jlo,i); j <= pert->jhi; j++) + fprintf(screen, "###FEP %2d %2d %9.5f %9.5f\n", i, j, + pert->array_orig[i][j], pert->array[i][j]); + } +#endif + + } else if (pert->which == ATOM) { + + if (pert->aparam == CHARGE) { // modify charges + int *atype = atom->type; + double *q = atom->q; + int *mask = atom->mask; + int natom = atom->nlocal + atom->nghost; + + for (i = 0; i < natom; i++) + if (atype[i] >= pert->ilo && atype[i] <= pert->ihi) + if (mask[i] & groupbit) + q[i] += delta; + +#ifdef FEP_MAXDEBUG + if (comm->me == 0 && screen) { + fprintf(screen, "###FEP change charge, delta = %f\n", delta); + fprintf(screen, "###FEP atom I old_q new_q\n"); + for (i = 0; i < atom->nlocal; i++) + if (atype[i] >= pert->ilo && atype[i] <= pert->ihi) + if (mask[i] & groupbit) + fprintf(screen, "###FEP %5d %2d %9.5f %9.5f\n", i, atype[i], + q_orig[i], q[i]); + } +#endif + + } + } + } + + // re-initialize pair styles if any PAIR settings were changed + // this resets other coeffs that may depend on changed values, + // and also offset and tail corrections + + if (pairflag) force->pair->reinit(); +} + + +/* ---------------------------------------------------------------------- + backup pair parameters +------------------------------------------------------------------------- */ + +void ComputeFEP::backup_params() +{ + int i,j; + + for (int m = 0; m < npert; m++) { + Perturb *pert = &perturb[m]; + if (pert->which == PAIR) { + for (i = pert->ilo; i <= pert->ihi; i++) + for (j = MAX(pert->jlo,i); j <= pert->jhi; j++) + pert->array_orig[i][j] = pert->array[i][j]; + } + } +} + + +/* ---------------------------------------------------------------------- + restore pair parameters to original values +------------------------------------------------------------------------- */ + +void ComputeFEP::restore_params() +{ + int i,j; + + for (int m = 0; m < npert; m++) { + Perturb *pert = &perturb[m]; + if (pert->which == PAIR) { + for (i = pert->ilo; i <= pert->ihi; i++) + for (j = MAX(pert->jlo,i); j <= pert->jhi; j++) + pert->array[i][j] = pert->array_orig[i][j]; + +#ifdef FEP_MAXDEBUG + if (comm->me == 0 && screen) { + fprintf(screen, "###FEP restore %s %s\n", pert->pstyle, pert->pparam); + fprintf(screen, "###FEP I J param\n"); + for (i = pert->ilo; i <= pert->ihi; i++) + for (j = MAX(pert->jlo,i); j <= pert->jhi; j++) + fprintf(screen, "###FEP %2d %2d %9.5f\n", i, j, pert->array[i][j]); + } +#endif + } +#ifdef FEP_MAXDEBUG + if (pert->which == ATOM) { + if (comm->me == 0 && screen) { + int *atype = atom->type; + double *q = atom->q; + int *mask = atom->mask; + int natom = atom->nlocal; + fprintf(screen, "###FEP restore charge\n"); + fprintf(screen, "###FEP atom I q\n"); + for (i = 0; i < natom; i++) + if (atype[i] >= pert->ilo && atype[i] <= pert->ihi) + if (mask[i] & groupbit) + fprintf(screen, "###FEP %5d %2d %9.5f\n", i, atype[i], q[i]); + } + } +#endif + } + + // re-initialize pair styles if any PAIR settings were changed + // this resets other coeffs that may depend on changed values, + // and also offset and tail corrections + + if (pairflag) force->pair->reinit(); +} + + +/* ---------------------------------------------------------------------- + manage storage for charge, force, energy, virial arrays +------------------------------------------------------------------------- */ + +void ComputeFEP::allocate_storage() +{ + nmax = atom->nmax; + if (chgflag) + memory->create(q_orig,nmax,"fep:q_orig"); + memory->create(f_orig,nmax,3,"fep:f_orig"); + memory->create(peatom_orig,nmax,"fep:peatom_orig"); + memory->create(pvatom_orig,nmax,6,"fep:pvatom_orig"); + if (force->kspace) { + memory->create(keatom_orig,nmax,"fep:keatom_orig"); + memory->create(kvatom_orig,nmax,6,"fep:kvatom_orig"); + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFEP::deallocate_storage() +{ + if (chgflag) + memory->destroy(q_orig); + memory->destroy(f_orig); + memory->destroy(peatom_orig); + memory->destroy(pvatom_orig); + if (force->kspace) { + memory->destroy(keatom_orig); + memory->destroy(kvatom_orig); + } +} + + +/* ---------------------------------------------------------------------- + backup and restore arrays with charge, force, energy, virial +------------------------------------------------------------------------- */ + +void ComputeFEP::backup_qfev() +{ + int i; + + int nall = atom->nlocal + atom->nghost; + int natom = atom->nlocal; + if (force->newton || force->kspace->tip4pflag) + natom += atom->nghost; + + if (chgflag) { + double *q = atom->q; + for (i = 0; i < nall; i++) + q_orig[i] = q[i]; + } + + double **f = atom->f; + for (i = 0; i < natom; i++) { + f_orig[i][0] = f[i][0]; + f_orig[i][1] = f[i][1]; + f_orig[i][2] = f[i][2]; + } + + eng_vdwl_orig = force->pair->eng_vdwl; + eng_coul_orig = force->pair->eng_coul; + + pvirial_orig[0] = force->pair->virial[0]; + pvirial_orig[1] = force->pair->virial[1]; + pvirial_orig[2] = force->pair->virial[2]; + pvirial_orig[3] = force->pair->virial[3]; + pvirial_orig[4] = force->pair->virial[4]; + pvirial_orig[5] = force->pair->virial[5]; + + if (update->eflag_atom) { + double *peatom = force->pair->eatom; + for (i = 0; i < natom; i++) + peatom_orig[i] = peatom[i]; + } + if (update->vflag_atom) { + double **pvatom = force->pair->vatom; + for (i = 0; i < natom; i++) { + pvatom_orig[i][0] = pvatom[i][0]; + pvatom_orig[i][1] = pvatom[i][1]; + pvatom_orig[i][2] = pvatom[i][2]; + pvatom_orig[i][3] = pvatom[i][3]; + pvatom_orig[i][4] = pvatom[i][4]; + pvatom_orig[i][5] = pvatom[i][5]; + } + } + + if (force->kspace) { + energy_orig = force->kspace->energy; + kvirial_orig[0] = force->kspace->virial[0]; + kvirial_orig[1] = force->kspace->virial[1]; + kvirial_orig[2] = force->kspace->virial[2]; + kvirial_orig[3] = force->kspace->virial[3]; + kvirial_orig[4] = force->kspace->virial[4]; + kvirial_orig[5] = force->kspace->virial[5]; + + if (update->eflag_atom) { + double *keatom = force->kspace->eatom; + for (i = 0; i < natom; i++) + keatom_orig[i] = keatom[i]; + } + if (update->vflag_atom) { + double **kvatom = force->kspace->vatom; + for (i = 0; i < natom; i++) { + kvatom_orig[i][0] = kvatom[i][0]; + kvatom_orig[i][1] = kvatom[i][1]; + kvatom_orig[i][2] = kvatom[i][2]; + kvatom_orig[i][3] = kvatom[i][3]; + kvatom_orig[i][4] = kvatom[i][4]; + kvatom_orig[i][5] = kvatom[i][5]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFEP::restore_qfev() +{ + int i; + + int nall = atom->nlocal + atom->nghost; + int natom = atom->nlocal; + if (force->newton || force->kspace->tip4pflag) + natom += atom->nghost; + + if (chgflag) { + double *q = atom->q; + for (i = 0; i < nall; i++) + q[i] = q_orig[i]; + } + + double **f = atom->f; + for (i = 0; i < natom; i++) { + f[i][0] = f_orig[i][0]; + f[i][1] = f_orig[i][1]; + f[i][2] = f_orig[i][2]; + } + + force->pair->eng_vdwl = eng_vdwl_orig; + force->pair->eng_coul = eng_coul_orig; + + force->pair->virial[0] = pvirial_orig[0]; + force->pair->virial[1] = pvirial_orig[1]; + force->pair->virial[2] = pvirial_orig[2]; + force->pair->virial[3] = pvirial_orig[3]; + force->pair->virial[4] = pvirial_orig[4]; + force->pair->virial[5] = pvirial_orig[5]; + + if (update->eflag_atom) { + double *peatom = force->pair->eatom; + for (i = 0; i < natom; i++) + peatom[i] = peatom_orig[i]; + } + if (update->vflag_atom) { + double **pvatom = force->pair->vatom; + for (i = 0; i < natom; i++) { + pvatom[i][0] = pvatom_orig[i][0]; + pvatom[i][1] = pvatom_orig[i][1]; + pvatom[i][2] = pvatom_orig[i][2]; + pvatom[i][3] = pvatom_orig[i][3]; + pvatom[i][4] = pvatom_orig[i][4]; + pvatom[i][5] = pvatom_orig[i][5]; + } + } + + if (force->kspace) { + force->kspace->energy = energy_orig; + force->kspace->virial[0] = kvirial_orig[0]; + force->kspace->virial[1] = kvirial_orig[1]; + force->kspace->virial[2] = kvirial_orig[2]; + force->kspace->virial[3] = kvirial_orig[3]; + force->kspace->virial[4] = kvirial_orig[4]; + force->kspace->virial[5] = kvirial_orig[5]; + + if (update->eflag_atom) { + double *keatom = force->kspace->eatom; + for (i = 0; i < natom; i++) + keatom[i] = keatom_orig[i]; + } + if (update->vflag_atom) { + double **kvatom = force->kspace->vatom; + for (i = 0; i < natom; i++) { + kvatom[i][0] = kvatom_orig[i][0]; + kvatom[i][1] = kvatom_orig[i][1]; + kvatom[i][2] = kvatom_orig[i][2]; + kvatom[i][3] = kvatom_orig[i][3]; + kvatom[i][4] = kvatom_orig[i][4]; + kvatom[i][5] = kvatom_orig[i][5]; + } + } + } +} + diff --git a/src/USER-FEP/compute_fep.h b/src/USER-FEP/compute_fep.h new file mode 100644 index 0000000000..7cc96034e7 --- /dev/null +++ b/src/USER-FEP/compute_fep.h @@ -0,0 +1,109 @@ +/* ---------------------------------------------------------------------- + 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: Agilio Padua (ICCF,UBP,CNRS) +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(fep,ComputeFEP) + +#else + +#ifndef COMPUTE_FEP_H +#define COMPUTE_FEP_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeFEP : public Compute { + public: + ComputeFEP(class LAMMPS *, int, char **); + ~ComputeFEP(); + void init(); + void compute_vector(); + + private: + int npert; + int pairflag; + int chgflag; + int tailflag, volumeflag; + int fepinitflag; + double temp_fep; + + int nmax; + double *q_orig; + double **f_orig; + double eng_vdwl_orig,eng_coul_orig; + double pvirial_orig[6]; + double *peatom_orig,**pvatom_orig; + double energy_orig; + double kvirial_orig[6]; + double *keatom_orig,**kvatom_orig; + + struct Perturb { + int which,ivar; + char *var; + char *pstyle,*pparam; + int ilo,ihi,jlo,jhi; + int pdim; + double **array,**array_orig; + int aparam; + }; + + Perturb *perturb; + + double compute_epair(); + void perturb_params(); + void backup_params(); + void restore_params(); + void allocate_storage(); + void deallocate_storage(); + void backup_qfev(); + void restore_qfev(); +}; + +} + +#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: Variable name for compute fep does not exist + +Self-explanatory. + +E: Variable for compute fep is invalid style + +Self-explanatory. + +E: Compute fep pair style does not exist + +Self-explanatory. + +E: Energy was not tallied on needed timestep + +You are using a thermo keyword that requires potentials to +have tallied energy, but they didn't on this timestep. See the +variable doc page for ideas on how to make this work. + +*/ diff --git a/src/USER-FEP/fix_adapt_fep.cpp b/src/USER-FEP/fix_adapt_fep.cpp new file mode 100644 index 0000000000..c2364889eb --- /dev/null +++ b/src/USER-FEP/fix_adapt_fep.cpp @@ -0,0 +1,460 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Charges by type and after option: Agilio Padua (Univ Blaise Pascal & CNRS) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "stdlib.h" +#include "fix_adapt_fep.h" +#include "atom.h" +#include "comm.h" +#include "update.h" +#include "modify.h" +#include "force.h" +#include "pair.h" +#include "pair_hybrid.h" +#include "kspace.h" +#include "input.h" +#include "variable.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace MathConst; + +enum{PAIR,KSPACE,ATOM}; +enum{DIAMETER,CHARGE}; + +#undef ADAPT_DEBUG + +/* ---------------------------------------------------------------------- */ + +FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 5) error->all(FLERR,"Illegal fix adapt/fep command"); + nevery = force->inumeric(FLERR,arg[3]); + if (nevery < 0) error->all(FLERR,"Illegal fix adapt/fep command"); + + // count # of adaptations + + nadapt = 0; + + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"pair") == 0) { + if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt/fep command"); + nadapt++; + iarg += 6; + } else if (strcmp(arg[iarg],"kspace") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command"); + nadapt++; + iarg += 2; + } else if (strcmp(arg[iarg],"atom") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix adapt/fep command"); + nadapt++; + iarg += 4; + } else break; + } + + if (nadapt == 0) error->all(FLERR,"Illegal fix adapt/fep command"); + adapt = new Adapt[nadapt]; + + // parse keywords + + nadapt = 0; + diamflag = 0; + chgflag = 0; + + iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"pair") == 0) { + if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt/fep command"); + adapt[nadapt].which = PAIR; + int n = strlen(arg[iarg+1]) + 1; + adapt[nadapt].pstyle = new char[n]; + strcpy(adapt[nadapt].pstyle,arg[iarg+1]); + n = strlen(arg[iarg+2]) + 1; + adapt[nadapt].pparam = new char[n]; + strcpy(adapt[nadapt].pparam,arg[iarg+2]); + force->bounds(arg[iarg+3],atom->ntypes, + adapt[nadapt].ilo,adapt[nadapt].ihi); + force->bounds(arg[iarg+4],atom->ntypes, + adapt[nadapt].jlo,adapt[nadapt].jhi); + if (strstr(arg[iarg+5],"v_") == arg[iarg+5]) { + n = strlen(&arg[iarg+5][2]) + 1; + adapt[nadapt].var = new char[n]; + strcpy(adapt[nadapt].var,&arg[iarg+5][2]); + } else error->all(FLERR,"Illegal fix adapt/fep command"); + nadapt++; + iarg += 6; + } else if (strcmp(arg[iarg],"kspace") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command"); + adapt[nadapt].which = KSPACE; + if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { + int n = strlen(&arg[iarg+1][2]) + 1; + adapt[nadapt].var = new char[n]; + strcpy(adapt[nadapt].var,&arg[iarg+1][2]); + } else error->all(FLERR,"Illegal fix adapt/fep command"); + nadapt++; + iarg += 2; + } else if (strcmp(arg[iarg],"atom") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt/fep command"); + adapt[nadapt].which = ATOM; + if (strcmp(arg[iarg+1],"diameter") == 0) { + adapt[nadapt].aparam = DIAMETER; + diamflag = 1; + } else if (strcmp(arg[iarg+1],"charge") == 0) { + adapt[nadapt].aparam = CHARGE; + chgflag = 1; + } else error->all(FLERR,"Illegal fix adapt/fep command"); + force->bounds(arg[iarg+2],atom->ntypes, + adapt[nadapt].ilo,adapt[nadapt].ihi); + if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { + int n = strlen(&arg[iarg+3][2]) + 1; + adapt[nadapt].var = new char[n]; + strcpy(adapt[nadapt].var,&arg[iarg+3][2]); + } else error->all(FLERR,"Illegal fix adapt/fep command"); + nadapt++; + iarg += 4; + } else break; + } + + // optional keywords + + resetflag = 0; + scaleflag = 0; + afterflag = 0; + + while (iarg < narg) { + if (strcmp(arg[iarg],"reset") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command"); + if (strcmp(arg[iarg+1],"no") == 0) resetflag = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) resetflag = 1; + else error->all(FLERR,"Illegal fix adapt/fep command"); + iarg += 2; + } else if (strcmp(arg[iarg],"scale") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command"); + if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1; + else error->all(FLERR,"Illegal fix adapt/fep command"); + iarg += 2; + } else if (strcmp(arg[iarg],"after") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command"); + if (strcmp(arg[iarg+1],"no") == 0) afterflag = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) afterflag = 1; + else error->all(FLERR,"Illegal fix adapt/fep command"); + iarg += 2; + } else error->all(FLERR,"Illegal fix adapt/fep command"); + } + + // allocate pair style arrays + + int n = atom->ntypes; + for (int m = 0; m < nadapt; m++) + if (adapt[m].which == PAIR) + memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig"); +} + +/* ---------------------------------------------------------------------- */ + +FixAdaptFEP::~FixAdaptFEP() +{ + for (int m = 0; m < nadapt; m++) { + delete [] adapt[m].var; + if (adapt[m].which == PAIR) { + delete [] adapt[m].pstyle; + delete [] adapt[m].pparam; + memory->destroy(adapt[m].array_orig); + } + } + delete [] adapt; +} + +/* ---------------------------------------------------------------------- */ + +int FixAdaptFEP::setmask() +{ + int mask = 0; + mask |= PRE_FORCE; + mask |= POST_RUN; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixAdaptFEP::init() +{ + int i,j; + + // setup and error checks + + anypair = 0; + + for (int m = 0; m < nadapt; m++) { + Adapt *ad = &adapt[m]; + + ad->ivar = input->variable->find(ad->var); + if (ad->ivar < 0) + error->all(FLERR,"Variable name for fix adapt/fep does not exist"); + if (!input->variable->equalstyle(ad->ivar)) + error->all(FLERR,"Variable for fix adapt/fep is invalid style"); + + if (ad->which == PAIR) { + anypair = 1; + + Pair *pair = force->pair_match(ad->pstyle,1); + if (pair == NULL) error->all(FLERR,"Fix adapt/fep pair style does not exist"); + void *ptr = pair->extract(ad->pparam,ad->pdim); + if (ptr == NULL) + error->all(FLERR,"Fix adapt/fep pair style param not supported"); + + ad->pdim = 2; + if (ad->pdim == 0) ad->scalar = (double *) ptr; + if (ad->pdim == 2) ad->array = (double **) ptr; + + // if pair hybrid, test that ilo,ihi,jlo,jhi are valid for sub-style + + if (ad->pdim == 2 && (strcmp(force->pair_style,"hybrid") == 0 || + strcmp(force->pair_style,"hybrid/overlay") == 0)) { + PairHybrid *pair = (PairHybrid *) force->pair; + for (i = ad->ilo; i <= ad->ihi; i++) + for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) + if (!pair->check_ijtype(i,j,ad->pstyle)) + error->all(FLERR,"Fix adapt/fep type pair range is not valid for " + "pair hybrid sub-style"); + } + + } else if (ad->which == KSPACE) { + if (force->kspace == NULL) + error->all(FLERR,"Fix adapt/fep kspace style does not exist"); + kspace_scale = (double *) force->kspace->extract("scale"); + + } else if (ad->which == ATOM) { + if (ad->aparam == DIAMETER) { + if (!atom->radius_flag) + error->all(FLERR,"Fix adapt/fep requires atom attribute diameter"); + } + if (ad->aparam == CHARGE) { + if (!atom->q_flag) + error->all(FLERR,"Fix adapt/fep requires atom attribute charge"); + } + } + } + + // make copy of original pair array values + + for (int m = 0; m < nadapt; m++) { + Adapt *ad = &adapt[m]; + if (ad->which == PAIR && ad->pdim == 2) { + for (i = ad->ilo; i <= ad->ihi; i++) + for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) + ad->array_orig[i][j] = ad->array[i][j]; + } + } + +#ifdef ADAPT_DEBUG + if (comm->me == 0 && screen) { + for (int m = 0; m < nadapt; m++) { + Adapt *ad = &adapt[m]; + if (ad->which == PAIR && ad->pdim == 2) { + fprintf(screen, "###ADAPT original %s %s\n", ad->pstyle, ad->pparam); + fprintf(screen, "###ADAPT I J old_param\n"); + for (i = ad->ilo; i <= ad->ihi; i++) + for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) + fprintf(screen, "###ADAPT %2d %2d %9.5f\n", i, j, + ad->array_orig[i][j]); + } + } + } +#endif + +} + +/* ---------------------------------------------------------------------- */ + +void FixAdaptFEP::setup_pre_force(int vflag) +{ + change_settings(); +} + +/* ---------------------------------------------------------------------- */ + +void FixAdaptFEP::pre_force(int vflag) +{ + if (nevery == 0) return; + + if (afterflag) { // update at n+1 (better with fix ave/time) + if (nevery == 1 || update->ntimestep == 0) + change_settings(); + else if (update->ntimestep > 1 && !((update->ntimestep - 1) % nevery)) + change_settings(); + } else { // original version: update at n + if (update->ntimestep % nevery) + return; + change_settings(); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixAdaptFEP::post_run() +{ + if (resetflag) restore_settings(); +} + +/* ---------------------------------------------------------------------- + change pair,kspace,atom parameters based on variable evaluation +------------------------------------------------------------------------- */ + +void FixAdaptFEP::change_settings() +{ + int i,j; + + // variable evaluation may invoke computes so wrap with clear/add + + modify->clearstep_compute(); + + for (int m = 0; m < nadapt; m++) { + Adapt *ad = &adapt[m]; + double value = input->variable->compute_equal(ad->ivar); + + // set global scalar or type pair array values + + if (ad->which == PAIR) { + if (ad->pdim == 0) { + if (scaleflag) *ad->scalar = value * ad->scalar_orig; + else *ad->scalar = value; + } else if (ad->pdim == 2) { + if (scaleflag) + for (i = ad->ilo; i <= ad->ihi; i++) + for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) + ad->array[i][j] = value*ad->array_orig[i][j]; + else + for (i = ad->ilo; i <= ad->ihi; i++) + for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) + ad->array[i][j] = value; + +#ifdef ADAPT_DEBUG + if (comm->me == 0 && screen) { + fprintf(screen, "###ADAPT change %s %s\n", ad->pstyle, ad->pparam); + fprintf(screen, "###ADAPT I J param\n"); + for (i = ad->ilo; i <= ad->ihi; i++) + for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) + fprintf(screen, "###ADAPT %2d %2d %9.5f\n", i, j, ad->array[i][j]); + } +#endif + + } + + // set kspace scale factor + + } else if (ad->which == KSPACE) { + *kspace_scale = value; + + } else if (ad->which == ATOM) { + + // set radius from diameter + // also scale rmass to new value + + if (ad->aparam == DIAMETER) { + int mflag = 0; + if (atom->rmass_flag) mflag = 1; + double density; + + int *atype = atom->type; + double *radius = atom->radius; + double *rmass = atom->rmass; + int *mask = atom->mask; + int natom = atom->nlocal + atom->nghost; + + if (mflag == 0) { + for (i = 0; i < natom; i++) + if (atype[i] >= ad->ilo && atype[i] <= ad->ihi) + if (mask[i] & groupbit) + radius[i] = 0.5*value; + } else { + for (i = 0; i < natom; i++) + if (atype[i] >= ad->ilo && atype[i] <= ad->ihi) + if (mask[i] & groupbit) { + density = rmass[i] / (4.0*MY_PI/3.0 * + radius[i]*radius[i]*radius[i]); + radius[i] = 0.5*value; + rmass[i] = 4.0*MY_PI/3.0 * + radius[i]*radius[i]*radius[i] * density; + } + } + } else if (ad->aparam == CHARGE) { + int *atype = atom->type; + double *q = atom->q; + int *mask = atom->mask; + int natom = atom->nlocal + atom->nghost; + + for (i = 0; i < natom; i++) + if (atype[i] >= ad->ilo && atype[i] <= ad->ihi) + if (mask[i] & groupbit) q[i] = value; + +#ifdef ADAPT_DEBUG + if (comm->me == 0 && screen) { + fprintf(screen, "###ADAPT change charge\n"); + fprintf(screen, "###ADAPT atom I q\n"); + for (i = 0; i < atom->nlocal; i++) + if (atype[i] >= ad->ilo && atype[i] <= ad->ihi) + if (mask[i] & groupbit) + fprintf(screen, "###ADAPT %5d %2d %9.5f\n", i, atype[i], q[i]); + } +#endif + } + } + } + + modify->addstep_compute(update->ntimestep + nevery); + + // re-initialize pair styles if any PAIR settings were changed + // this resets other coeffs that may depend on changed values, + // and also offset and tail corrections + + if (anypair) force->pair->reinit(); +} + +/* ---------------------------------------------------------------------- + restore pair,kspace.atom parameters to original values +------------------------------------------------------------------------- */ + +void FixAdaptFEP::restore_settings() +{ + for (int m = 0; m < nadapt; m++) { + Adapt *ad = &adapt[m]; + if (ad->which == PAIR) { + if (ad->pdim == 0) *ad->scalar = ad->scalar_orig; + else if (ad->pdim == 2) { + for (int i = ad->ilo; i <= ad->ihi; i++) + for (int j = MAX(ad->jlo,i); j <= ad->jhi; j++) + ad->array[i][j] = ad->array_orig[i][j]; + } + + } else if (ad->which == KSPACE) { + *kspace_scale = 1.0; + + } else if (ad->which == ATOM) { + + } + } + + if (anypair) force->pair->reinit(); +} diff --git a/src/USER-FEP/fix_adapt_fep.h b/src/USER-FEP/fix_adapt_fep.h new file mode 100644 index 0000000000..0fca92fcaa --- /dev/null +++ b/src/USER-FEP/fix_adapt_fep.h @@ -0,0 +1,107 @@ +/* ---------------------------------------------------------------------- + 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(adapt/fep,FixAdaptFEP) + +#else + +#ifndef LMP_FIX_ADAPT_FEP_H +#define LMP_FIX_ADAPT_FEP_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixAdaptFEP : public Fix { + public: + int diamflag; // 1 if atom diameters will vary, for AtomVecGranular + int chgflag; + + FixAdaptFEP(class LAMMPS *, int, char **); + ~FixAdaptFEP(); + int setmask(); + void init(); + void setup_pre_force(int); + void pre_force(int); + void post_run(); + + private: + int nadapt,resetflag,scaleflag,afterflag; + int anypair; + + struct Adapt { + int which,ivar; + char *var; + char *pstyle,*pparam; + int ilo,ihi,jlo,jhi; + int pdim; + double *scalar,scalar_orig; + double **array,**array_orig; + int aparam; + }; + + Adapt *adapt; + double *kspace_scale; + + void change_settings(); + void restore_settings(); +}; + +} + +#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: Variable name for fix adapt does not exist + +Self-explanatory. + +E: Variable for fix adapt is invalid style + +Only equal-style variables can be used. + +E: Fix adapt pair style does not exist + +Self-explanatory + +E: Fix adapt pair style param not supported + +The pair style does not know about the parameter you specified. + +E: Fix adapt type pair range is not valid for pair hybrid sub-style + +Self-explanatory. + +E: Fix adapt kspace style does not exist + +Self-explanatory. + +E: Fix adapt requires atom attribute diameter + +The atom style being used does not specify an atom diameter. + +E: Fix adapt requires atom attribute charge + +The atom style being used does not specify an atom charge. + +*/ diff --git a/src/USER-FEP/pair_coul_cut_soft.cpp b/src/USER-FEP/pair_coul_cut_soft.cpp new file mode 100644 index 0000000000..831f212c70 --- /dev/null +++ b/src/USER-FEP/pair_coul_cut_soft.cpp @@ -0,0 +1,376 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_coul_cut_soft.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairCoulCutSoft::PairCoulCutSoft(LAMMPS *lmp) : Pair(lmp) {} + +/* ---------------------------------------------------------------------- */ + +PairCoulCutSoft::~PairCoulCutSoft() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(lambda); + memory->destroy(lam1); + memory->destroy(lam2); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairCoulCutSoft::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; + double rsq,forcecoul,factor_coul; + double denc; + int *ilist,*jlist,*numneigh,**firstneigh; + + 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; + + 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]; + qtmp = q[i]; + 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; + + 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]) { + + denc = sqrt(lam2[itype][jtype] + rsq); + forcecoul = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / (denc*denc*denc); + + fpair = factor_coul*forcecoul; + + 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_coul * qqrd2e * lam1[itype][jtype] * qtmp*q[j] / denc; + + 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 PairCoulCutSoft::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(lambda,n+1,n+1,"pair:lambda"); + memory->create(lam1,n+1,n+1,"pair:lam1"); + memory->create(lam2,n+1,n+1,"pair:lam2"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairCoulCutSoft::settings(int narg, char **arg) +{ + if (narg != 3) error->all(FLERR,"Illegal pair_style command"); + + nlambda = force->numeric(FLERR,arg[0]); + alphac = force->numeric(FLERR,arg[1]); + + cut_global = 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]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairCoulCutSoft::coeff(int narg, char **arg) +{ + if (narg < 3 || narg > 4) + 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 lambda_one = force->numeric(FLERR,arg[2]); + + double cut_one = cut_global; + if (narg == 4) cut_one = force->numeric(FLERR,arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + lambda[i][j] = lambda_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairCoulCutSoft::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style coul/cut/soft requires atom attribute q"); + + neighbor->request(this); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairCoulCutSoft::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + if (lambda[i][i] != lambda[j][j]) + error->all(FLERR,"Pair coul/cut/soft different lambda values in mix"); + lambda[i][j] = lambda[i][i]; + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + lam1[i][j] = pow(lambda[i][j], nlambda); + lam2[i][j] = alphac * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); + + cut[j][i] = cut[i][j]; + lambda[j][i] = lambda[i][j]; + lam1[j][i] = lam1[i][j]; + lam2[j][i] = lam2[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairCoulCutSoft::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(&lambda[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairCoulCutSoft::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(&lambda[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&lambda[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairCoulCutSoft::write_restart_settings(FILE *fp) +{ + fwrite(&nlambda,sizeof(double),1,fp); + fwrite(&alphac,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 PairCoulCutSoft::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&nlambda,sizeof(double),1,fp); + fread(&alphac,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(&nlambda,1,MPI_DOUBLE,0,world); + MPI_Bcast(&alphac,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); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairCoulCutSoft::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g\n",i,lambda[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairCoulCutSoft::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\n",i,j,lambda[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairCoulCutSoft::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, + double &fforce) +{ + double forcecoul,phicoul; + double denc; + + if (rsq < cutsq[itype][jtype]) { + denc = sqrt(lam2[itype][jtype] + rsq); + forcecoul = force->qqrd2e * lam1[itype][jtype] * atom->q[i]*atom->q[j] / + (denc*denc*denc); + } else forcecoul = 0.0; + fforce = factor_coul*forcecoul; + + if (rsq < cutsq[itype][jtype]) + phicoul = force->qqrd2e * lam1[itype][jtype] * atom->q[i]*atom->q[j] / denc; + else phicoul = 0.0; + return factor_coul*phicoul; +} + +/* ---------------------------------------------------------------------- */ + +void *PairCoulCutSoft::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"lambda") == 0) return (void *) lambda; + return NULL; +} diff --git a/src/USER-FEP/pair_coul_cut_soft.h b/src/USER-FEP/pair_coul_cut_soft.h new file mode 100644 index 0000000000..25c6bd4142 --- /dev/null +++ b/src/USER-FEP/pair_coul_cut_soft.h @@ -0,0 +1,80 @@ +/* ---------------------------------------------------------------------- + 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(coul/cut/soft,PairCoulCutSoft) + +#else + +#ifndef LMP_PAIR_COUL_CUT_SOFT_H +#define LMP_PAIR_COUL_CUT_SOFT_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairCoulCutSoft : public Pair { + public: + PairCoulCutSoft(class LAMMPS *); + virtual ~PairCoulCutSoft(); + 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 *); + void write_data(FILE *); + void write_data_all(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); + void *extract(const char *, int &); + + protected: + double cut_global; + double **cut; + double **lambda; + double nlambda, alphac; + double **lam1, **lam2; + + 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 coul/cut/soft requires atom attribute q + +The atom style defined does not have these attributes. + +E: Pair coul/cut/soft different lambda values in mix + +The value of lambda has to be the same for I J pairs. + +*/ diff --git a/src/USER-FEP/pair_coul_long_soft.cpp b/src/USER-FEP/pair_coul_long_soft.cpp new file mode 100644 index 0000000000..d396744005 --- /dev/null +++ b/src/USER-FEP/pair_coul_long_soft.cpp @@ -0,0 +1,395 @@ +/* ---------------------------------------------------------------------- + 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) + Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_coul_long_soft.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "update.h" +#include "integrate.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairCoulLongSoft::PairCoulLongSoft(LAMMPS *lmp) : Pair(lmp) +{ + ewaldflag = pppmflag = 1; + qdist = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +PairCoulLongSoft::~PairCoulLongSoft() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(scale); + + memory->destroy(lambda); + memory->destroy(lam1); + memory->destroy(lam2); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairCoulLongSoft::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; + double r,rsq,forcecoul,factor_coul; + double grij,expm2,prefactor,t,erfc; + double denc; + int *ilist,*jlist,*numneigh,**firstneigh; + + 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; + + 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]; + qtmp = q[i]; + 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; + + 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 < cut_coulsq) { + + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + + denc = sqrt(lam2[itype][jtype] + rsq); + prefactor = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / (denc*denc*denc); + + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + + fpair = forcecoul; + + 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) { + prefactor = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / denc; + ecoul = prefactor*erfc; + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } + + 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 PairCoulLongSoft::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(scale,n+1,n+1,"pair:scale"); + + memory->create(lambda,n+1,n+1,"pair:lambda"); + memory->create(lam1,n+1,n+1,"pair:lam1"); + memory->create(lam2,n+1,n+1,"pair:lam2"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairCoulLongSoft::settings(int narg, char **arg) +{ + if (narg != 3) error->all(FLERR,"Illegal pair_style command"); + + nlambda = force->numeric(FLERR,arg[0]); + alphac = force->numeric(FLERR,arg[1]); + + cut_coul = force->numeric(FLERR,arg[2]); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairCoulLongSoft::coeff(int narg, char **arg) +{ + if (narg != 3) + 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 lambda_one = force->numeric(FLERR,arg[2]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + lambda[i][j] = lambda_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 PairCoulLongSoft::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style lj/cut/coul/long requires atom attribute q"); + + neighbor->request(this); + + cut_coulsq = cut_coul * cut_coul; + + // 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; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairCoulLongSoft::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + if (lambda[i][i] != lambda[j][j]) + error->all(FLERR,"Pair coul/cut/soft different lambda values in mix"); + lambda[i][j] = lambda[i][i]; + } + + lam1[i][j] = pow(lambda[i][j], nlambda); + lam2[i][j] = alphac * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); + + scale[j][i] = scale[i][j]; + lambda[j][i] = lambda[i][j]; + lam1[j][i] = lam1[i][j]; + lam2[j][i] = lam2[i][j]; + + return cut_coul+2.0*qdist; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairCoulLongSoft::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(&lambda[i][j],sizeof(double),1,fp); + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairCoulLongSoft::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(&lambda[i][j],sizeof(double),1,fp); + MPI_Bcast(&lambda[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairCoulLongSoft::write_restart_settings(FILE *fp) +{ + fwrite(&nlambda,sizeof(double),1,fp); + fwrite(&alphac,sizeof(double),1,fp); + + fwrite(&cut_coul,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 PairCoulLongSoft::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&nlambda,sizeof(double),1,fp); + fread(&alphac,sizeof(double),1,fp); + + fread(&cut_coul,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world); + MPI_Bcast(&alphac,1,MPI_DOUBLE,0,world); + + MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +double PairCoulLongSoft::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r,grij,expm2,t,erfc,prefactor; + double forcecoul,phicoul; + double denc; + + if (rsq < cut_coulsq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + + denc = sqrt(lam2[itype][jtype] + rsq); + prefactor = force->qqrd2e * lam1[itype][jtype] * atom->q[i]*atom->q[j] / + (denc*denc*denc); + + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + + fforce = forcecoul; + + if (rsq < cut_coulsq) { + prefactor = force->qqrd2e * lam1[itype][jtype] * atom->q[i]*atom->q[j] / denc; + phicoul = prefactor*erfc; + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + } else phicoul = 0.0; + + return phicoul; +} + +/* ---------------------------------------------------------------------- */ + +void *PairCoulLongSoft::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + dim = 2; + if (strcmp(str,"scale") == 0) return (void *) scale; + if (strcmp(str,"lambda") == 0) return (void *) lambda; + + return NULL; +} diff --git a/src/USER-FEP/pair_coul_long_soft.h b/src/USER-FEP/pair_coul_long_soft.h new file mode 100644 index 0000000000..013a94a2f4 --- /dev/null +++ b/src/USER-FEP/pair_coul_long_soft.h @@ -0,0 +1,80 @@ +/* -*- 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(coul/long/soft,PairCoulLongSoft) + +#else + +#ifndef LMP_PAIR_COUL_LONG_SOFT_H +#define LMP_PAIR_COUL_LONG_SOFT_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairCoulLongSoft : public Pair { + public: + PairCoulLongSoft(class LAMMPS *); + virtual ~PairCoulLongSoft(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + virtual 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 &); + virtual void *extract(const char *, int &); + + protected: + double cut_coul,cut_coulsq; + double **scale; + double **lambda; + double nlambda, alphac; + double **lam1, **lam2; + double qdist; // TIP4P distance O to negative charge (compatibility of cutoffs) + double g_ewald; + + 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. + +*/ diff --git a/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp b/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp new file mode 100644 index 0000000000..ac2e4b4af5 --- /dev/null +++ b/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp @@ -0,0 +1,514 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lj_cut_coul_cut_soft.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairLJCutCoulCutSoft::PairLJCutCoulCutSoft(LAMMPS *lmp) : Pair(lmp) +{ + writedata = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairLJCutCoulCutSoft::~PairLJCutCoulCutSoft() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut_lj); + memory->destroy(cut_ljsq); + memory->destroy(cut_coul); + memory->destroy(cut_coulsq); + memory->destroy(epsilon); + memory->destroy(sigma); + memory->destroy(lambda); + memory->destroy(lj1); + memory->destroy(lj2); + memory->destroy(lj3); + memory->destroy(lj4); + memory->destroy(offset); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulCutSoft::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double rsq,forcecoul,forcelj,factor_coul,factor_lj; + double denc, denlj, r4sig6; + int *ilist,*jlist,*numneigh,**firstneigh; + + 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; + + 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]; + qtmp = q[i]; + 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_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]) { + + if (rsq < cut_coulsq[itype][jtype]) { + denc = sqrt(lj4[itype][jtype] + rsq); + forcecoul = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc); + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + } else forcelj = 0.0; + + fpair = factor_coul*forcecoul + factor_lj*forcelj; + + 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[itype][jtype]) + ecoul = factor_coul * qqrd2e * lj1[itype][jtype] * qtmp*q[j] / denc; + else ecoul = 0.0; + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * + (1.0/(denlj*denlj) - 1.0/denlj) - 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 PairLJCutCoulCutSoft::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(cut_coul,n+1,n+1,"pair:cut_coul"); + memory->create(cut_coulsq,n+1,n+1,"pair:cut_coulsq"); + memory->create(epsilon,n+1,n+1,"pair:epsilon"); + memory->create(sigma,n+1,n+1,"pair:sigma"); + memory->create(lambda,n+1,n+1,"pair:lambda"); + 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 PairLJCutCoulCutSoft::settings(int narg, char **arg) +{ + if (narg < 4 || narg > 5) error->all(FLERR,"Illegal pair_style command"); + + nlambda = force->numeric(FLERR,arg[0]); + alphalj = force->numeric(FLERR,arg[1]); + alphac = force->numeric(FLERR,arg[2]); + + cut_lj_global = force->numeric(FLERR,arg[3]); + if (narg == 4) cut_coul_global = cut_lj_global; + else cut_coul_global = force->numeric(FLERR,arg[4]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) { + cut_lj[i][j] = cut_lj_global; + cut_coul[i][j] = cut_coul_global; + } + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJCutCoulCutSoft::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 lambda_one = force->numeric(FLERR,arg[4]); + + double cut_lj_one = cut_lj_global; + double cut_coul_one = cut_coul_global; + if (narg >= 6) cut_coul_one = cut_lj_one = force->numeric(FLERR,arg[5]); + if (narg == 7) cut_coul_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; + lambda[i][j] = lambda_one; + cut_lj[i][j] = cut_lj_one; + cut_coul[i][j] = cut_coul_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJCutCoulCutSoft::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style lj/cut/coul/cut/soft requires atom attribute q"); + + neighbor->request(this); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLJCutCoulCutSoft::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]); + if (lambda[i][i] != lambda[j][j]) + error->all(FLERR,"Pair lj/cut/coul/cut/soft different lambda values in mix"); + lambda[i][j] = lambda[i][i]; + cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); + cut_coul[i][j] = mix_distance(cut_coul[i][i],cut_coul[j][j]); + } + + double cut = MAX(cut_lj[i][j],cut_coul[i][j]); + cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; + cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j]; + + lj1[i][j] = pow(lambda[i][j], nlambda); + lj2[i][j] = pow(sigma[i][j], 6.0); + lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); + lj4[i][j] = alphac * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); + + if (offset_flag) { + double denlj = lj3[i][j] + pow(sigma[i][j] / cut_lj[i][j], 6.0); + offset[i][j] = lj1[i][j] * 4.0 * epsilon[i][j] * (1.0/(denlj*denlj) - 1.0/denlj); + } else offset[i][j] = 0.0; + + epsilon[j][i] = epsilon[i][j]; + sigma[j][i] = sigma[i][j]; + lambda[j][i] = lambda[i][j]; + cut_ljsq[j][i] = cut_ljsq[i][j]; + cut_coulsq[j][i] = cut_coulsq[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]; + + // 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]* lj1[i][j] * epsilon[i][j] * + sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); + ptail_ij = 16.0*MY_PI*all[0]*all[1]* lj1[i][j] * epsilon[i][j] * + sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); + } + + return cut; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJCutCoulCutSoft::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(&lambda[i][j],sizeof(double),1,fp); + fwrite(&cut_lj[i][j],sizeof(double),1,fp); + fwrite(&cut_coul[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutCoulCutSoft::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(&lambda[i][j],sizeof(double),1,fp); + fread(&cut_lj[i][j],sizeof(double),1,fp); + fread(&cut_coul[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(&lambda[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJCutCoulCutSoft::write_restart_settings(FILE *fp) +{ + fwrite(&nlambda,sizeof(double),1,fp); + fwrite(&alphalj,sizeof(double),1,fp); + fwrite(&alphac,sizeof(double),1,fp); + + fwrite(&cut_lj_global,sizeof(double),1,fp); + fwrite(&cut_coul_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); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutCoulCutSoft::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&nlambda,sizeof(double),1,fp); + fread(&alphalj,sizeof(double),1,fp); + fread(&alphac,sizeof(double),1,fp); + + fread(&cut_lj_global,sizeof(double),1,fp); + fread(&cut_coul_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); + } + + MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world); + MPI_Bcast(&alphalj,1,MPI_DOUBLE,0,world); + MPI_Bcast(&alphac,1,MPI_DOUBLE,0,world); + + MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul_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); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairLJCutCoulCutSoft::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],lambda[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairLJCutCoulCutSoft::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 %g\n",i,j,epsilon[i][j],sigma[i][j], + lambda[i][j],cut_lj[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJCutCoulCutSoft::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double forcecoul,forcelj,phicoul,philj; + double denc, denlj, r4sig6; + + if (rsq < cut_coulsq[itype][jtype]) { + denc = sqrt(lj4[itype][jtype] + rsq); + forcecoul = force->qqrd2e * lj1[itype][jtype] * atom->q[i]*atom->q[j] / + (denc*denc*denc); + } else forcecoul = 0.0; + if (rsq < cut_ljsq[itype][jtype]) { + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + } else forcelj = 0.0; + fforce = factor_coul*forcecoul + factor_lj*forcelj; + + double eng = 0.0; + if (rsq < cut_coulsq[itype][jtype]) { + phicoul = force->qqrd2e * lj1[itype][jtype] * atom->q[i]*atom->q[j] / denc; + eng += factor_coul*phicoul; + } + if (rsq < cut_ljsq[itype][jtype]) { + philj = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * + (1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype]; + eng += factor_lj*philj; + } + + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJCutCoulCutSoft::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str,"lambda") == 0) return (void *) lambda; + return NULL; +} diff --git a/src/USER-FEP/pair_lj_cut_coul_cut_soft.h b/src/USER-FEP/pair_lj_cut_coul_cut_soft.h new file mode 100644 index 0000000000..fa9e0ab71f --- /dev/null +++ b/src/USER-FEP/pair_lj_cut_coul_cut_soft.h @@ -0,0 +1,82 @@ +/* -*- 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/coul/cut/soft,PairLJCutCoulCutSoft) + +#else + +#ifndef LMP_PAIR_LJ_CUT_COUL_CUT_SOFT_H +#define LMP_PAIR_LJ_CUT_COUL_CUT_SOFT_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLJCutCoulCutSoft : public Pair { + public: + PairLJCutCoulCutSoft(class LAMMPS *); + virtual ~PairLJCutCoulCutSoft(); + 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 *); + void write_data(FILE *); + void write_data_all(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); + void *extract(const char *, int &); + + protected: + double cut_lj_global,cut_coul_global; + double **cut_lj,**cut_ljsq; + double **cut_coul,**cut_coulsq; + double **epsilon,**sigma, **lambda; + double nlambda, alphalj, alphac; + double **lj1,**lj2,**lj3,**lj4,**offset; + + 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/cut/soft requires atom attribute q + +The atom style defined does not have this attribute. + +E: Pair lj/cut/coul/cut/soft different lambda values in mix + +The value of lambda has to be the same for I J pairs. + + +*/ diff --git a/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp b/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp new file mode 100644 index 0000000000..8e599f2886 --- /dev/null +++ b/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp @@ -0,0 +1,955 @@ +/* ---------------------------------------------------------------------- + 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) + Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lj_cut_coul_long_soft.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "update.h" +#include "integrate.h" +#include "respa.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 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairLJCutCoulLongSoft::PairLJCutCoulLongSoft(LAMMPS *lmp) : Pair(lmp) +{ + ewaldflag = pppmflag = 1; + respa_enable = 1; + writedata = 1; + qdist = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +PairLJCutCoulLongSoft::~PairLJCutCoulLongSoft() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut_lj); + memory->destroy(cut_ljsq); + memory->destroy(epsilon); + memory->destroy(sigma); + memory->destroy(lambda); + memory->destroy(lj1); + memory->destroy(lj2); + memory->destroy(lj3); + memory->destroy(lj4); + memory->destroy(offset); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::compute(int eflag, int vflag) +{ + int i,ii,j,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double r,rsq,forcecoul,forcelj,factor_coul,factor_lj; + double grij,expm2,prefactor,t,erfc; + double denc, denlj, r4sig6; + int *ilist,*jlist,*numneigh,**firstneigh; + + 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; + + 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]; + qtmp = q[i]; + 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_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]) { + + if (rsq < cut_coulsq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + + denc = sqrt(lj4[itype][jtype] + rsq); + prefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc); + + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + } else forcelj = 0.0; + + fpair = forcecoul + factor_lj*forcelj; + + 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) { + prefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / denc; + ecoul = prefactor*erfc; + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * + (1.0/(denlj*denlj) - 1.0/denlj) - 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(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::compute_inner() +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; + double rsq,forcecoul,forcelj,factor_coul,factor_lj; + double rsw; + double denc, denlj, r4sig6; + int *ilist,*jlist,*numneigh,**firstneigh; + + 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; + + inum = listinner->inum; + ilist = listinner->ilist; + numneigh = listinner->numneigh; + firstneigh = listinner->firstneigh; + + double cut_out_on = cut_respa[0]; + double cut_out_off = cut_respa[1]; + + double cut_out_diff = cut_out_off - cut_out_on; + double cut_out_on_sq = cut_out_on*cut_out_on; + double cut_out_off_sq = cut_out_off*cut_out_off; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + 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_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; + + if (rsq < cut_out_off_sq) { + jtype = type[j]; + + denc = sqrt(lj4[itype][jtype] + rsq); + forcecoul = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; + + if (rsq < cut_ljsq[itype][jtype]) { + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + } else forcelj = 0.0; + + fpair = forcecoul + factor_lj*forcelj; + + if (rsq > cut_out_on_sq) { + rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; + fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); + } + + 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; + } + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::compute_middle() +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; + double rsq,forcecoul,forcelj,factor_coul,factor_lj; + double rsw; + double denc, denlj, r4sig6; + int *ilist,*jlist,*numneigh,**firstneigh; + + 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; + + inum = listmiddle->inum; + ilist = listmiddle->ilist; + numneigh = listmiddle->numneigh; + firstneigh = listmiddle->firstneigh; + + double cut_in_off = cut_respa[0]; + double cut_in_on = cut_respa[1]; + double cut_out_on = cut_respa[2]; + double cut_out_off = cut_respa[3]; + + double cut_in_diff = cut_in_on - cut_in_off; + double cut_out_diff = cut_out_off - cut_out_on; + double cut_in_off_sq = cut_in_off*cut_in_off; + double cut_in_on_sq = cut_in_on*cut_in_on; + double cut_out_on_sq = cut_out_on*cut_out_on; + double cut_out_off_sq = cut_out_off*cut_out_off; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + 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_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; + + if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { + jtype = type[j]; + + denc = sqrt(lj4[itype][jtype] + rsq); + forcecoul = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; + + if (rsq < cut_ljsq[itype][jtype]) { + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + } else forcelj = 0.0; + + fpair = forcecoul + factor_lj*forcelj; + + if (rsq < cut_in_on_sq) { + rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; + fpair *= rsw*rsw*(3.0 - 2.0*rsw); + } + if (rsq > cut_out_on_sq) { + rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; + fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); + } + + 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; + } + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::compute_outer(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double r,rsq,forcecoul,forcelj,factor_coul,factor_lj; + double grij,expm2,fprefactor,eprefactor,t,erfc; + double rsw; + double denc, denlj, r4sig6; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 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; + + inum = listouter->inum; + ilist = listouter->ilist; + numneigh = listouter->numneigh; + firstneigh = listouter->firstneigh; + + double cut_in_off = cut_respa[2]; + double cut_in_on = cut_respa[3]; + + double cut_in_diff = cut_in_on - cut_in_off; + double cut_in_off_sq = cut_in_off*cut_in_off; + double cut_in_on_sq = cut_in_on*cut_in_on; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + 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_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]) { + + if (rsq < cut_coulsq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + + denc = sqrt(lj4[itype][jtype] + rsq); + fprefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc); + + forcecoul = fprefactor * (erfc + EWALD_F*grij*expm2 - 1.0); + + if (rsq > cut_in_off_sq) { + if (rsq < cut_in_on_sq) { + rsw = (r - cut_in_off)/cut_in_diff; + forcecoul += fprefactor*rsw*rsw*(3.0 - 2.0*rsw); + if (factor_coul < 1.0) + forcecoul -= + (1.0-factor_coul)*fprefactor*rsw*rsw*(3.0 - 2.0*rsw); + } else { + forcecoul += fprefactor; + if (factor_coul < 1.0) + forcecoul -= (1.0-factor_coul)*fprefactor; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) { + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + if (rsq < cut_in_on_sq) { + rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; + forcelj *= rsw*rsw*(3.0 - 2.0*rsw); + } + } else forcelj = 0.0; + + fpair = forcecoul + forcelj; + + 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) { + eprefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / denc; + ecoul = eprefactor*erfc; + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*eprefactor; + } else ecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * + (1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (vflag) { + if (rsq < cut_coulsq) { + forcecoul = fprefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*fprefactor; + } else forcecoul = 0.0; + + if (rsq <= cut_in_off_sq) { + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + } else if (rsq < cut_in_on_sq) { + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + } + fpair = forcecoul + factor_lj*forcelj; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::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(epsilon,n+1,n+1,"pair:epsilon"); + memory->create(sigma,n+1,n+1,"pair:sigma"); + memory->create(lambda,n+1,n+1,"pair:lambda"); + 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 PairLJCutCoulLongSoft::settings(int narg, char **arg) +{ + if (narg < 4 || narg > 5) error->all(FLERR,"Illegal pair_style command"); + + nlambda = force->numeric(FLERR,arg[0]); + alphalj = force->numeric(FLERR,arg[1]); + alphac = force->numeric(FLERR,arg[2]); + + cut_lj_global = force->numeric(FLERR,arg[3]); + if (narg == 4) cut_coul = cut_lj_global; + else cut_coul = force->numeric(FLERR,arg[4]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::coeff(int narg, char **arg) +{ + if (narg < 5 || narg > 6) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double epsilon_one = force->numeric(FLERR,arg[2]); + double sigma_one = force->numeric(FLERR,arg[3]); + double lambda_one = force->numeric(FLERR,arg[4]); + + double cut_lj_one = cut_lj_global; + if (narg == 6) cut_lj_one = force->numeric(FLERR,arg[5]); + + 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; + lambda[i][j] = lambda_one; + cut_lj[i][j] = cut_lj_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style lj/cut/coul/long/soft requires atom attribute q"); + + // request regular or rRESPA neighbor lists + + int irequest; + + if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { + int respa = 0; + if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; + if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + + if (respa == 0) irequest = neighbor->request(this); + else if (respa == 1) { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } else { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 2; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respamiddle = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } + + } else irequest = neighbor->request(this); + + cut_coulsq = cut_coul * cut_coul; + + // set rRESPA cutoffs + + if (strstr(update->integrate_style,"respa") && + ((Respa *) update->integrate)->level_inner >= 0) + cut_respa = ((Respa *) update->integrate)->cutoff; + else 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; +} + +/* ---------------------------------------------------------------------- + neighbor callback to inform pair style of neighbor list to use + regular or rRESPA +------------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::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 PairLJCutCoulLongSoft::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]); + if (lambda[i][i] != lambda[j][j]) + error->all(FLERR,"Pair lj/cut/coul/long/soft different lambda values in mix"); + lambda[i][j] = lambda[i][i]; + cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); + } + + // 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] = pow(lambda[i][j], nlambda); + lj2[i][j] = pow(sigma[i][j], 6.0); + lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); + lj4[i][j] = alphac * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); + + if (offset_flag) { + double denlj = lj3[i][j] + pow(sigma[i][j] / cut_lj[i][j], 6.0); + offset[i][j] = lj1[i][j] * 4.0 * epsilon[i][j] * (1.0/(denlj*denlj) - 1.0/denlj); + } else offset[i][j] = 0.0; + + epsilon[j][i] = epsilon[i][j]; + sigma[j][i] = sigma[i][j]; + lambda[j][i] = lambda[i][j]; + 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]; + + // 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 PairLJCutCoulLongSoft::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(&lambda[i][j],sizeof(double),1,fp); + fwrite(&cut_lj[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::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(&lambda[i][j],sizeof(double),1,fp); + 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(&lambda[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 PairLJCutCoulLongSoft::write_restart_settings(FILE *fp) +{ + fwrite(&nlambda,sizeof(double),1,fp); + fwrite(&alphalj,sizeof(double),1,fp); + fwrite(&alphac,sizeof(double),1,fp); + + fwrite(&cut_lj_global,sizeof(double),1,fp); + fwrite(&cut_coul,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); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&nlambda,sizeof(double),1,fp); + fread(&alphalj,sizeof(double),1,fp); + fread(&alphac,sizeof(double),1,fp); + + fread(&cut_lj_global,sizeof(double),1,fp); + fread(&cut_coul,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); + } + + MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world); + MPI_Bcast(&alphalj,1,MPI_DOUBLE,0,world); + MPI_Bcast(&alphac,1,MPI_DOUBLE,0,world); + + MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul,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); +} + + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],lambda[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairLJCutCoulLongSoft::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 %g\n",i,j,epsilon[i][j],sigma[i][j], + lambda[i][j],cut_lj[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJCutCoulLongSoft::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r,grij,expm2,t,erfc,prefactor; + double forcecoul,forcelj,phicoul,philj; + double denc, denlj, r4sig6; + + if (rsq < cut_coulsq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + + denc = sqrt(lj4[itype][jtype] + rsq); + prefactor = force->qqrd2e * lj1[itype][jtype] * atom->q[i]*atom->q[j] / + (denc*denc*denc); + + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + } else forcelj = 0.0; + + fforce = forcecoul + factor_lj*forcelj; + + double eng = 0.0; + if (rsq < cut_coulsq) { + prefactor = force->qqrd2e * lj1[itype][jtype] * atom->q[i]*atom->q[j] / denc; + phicoul = prefactor*erfc; + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + eng += phicoul; + } + + if (rsq < cut_ljsq[itype][jtype]) { + philj = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * + (1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype]; + eng += factor_lj*philj; + } + + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJCutCoulLongSoft::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str,"lambda") == 0) return (void *) lambda; + return NULL; +} diff --git a/src/USER-FEP/pair_lj_cut_coul_long_soft.h b/src/USER-FEP/pair_lj_cut_coul_long_soft.h new file mode 100644 index 0000000000..a03be3814a --- /dev/null +++ b/src/USER-FEP/pair_lj_cut_coul_long_soft.h @@ -0,0 +1,94 @@ +/* -*- 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/coul/long/soft,PairLJCutCoulLongSoft) + +#else + +#ifndef LMP_PAIR_LJ_CUT_COUL_LONG_SOFT_H +#define LMP_PAIR_LJ_CUT_COUL_LONG_SOFT_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLJCutCoulLongSoft : public Pair { + public: + PairLJCutCoulLongSoft(class LAMMPS *); + virtual ~PairLJCutCoulLongSoft(); + 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 &); + + void compute_inner(); + void compute_middle(); + virtual void compute_outer(int, int); + + protected: + double cut_lj_global; + double **cut_lj,**cut_ljsq; + double cut_coul,cut_coulsq; + double **epsilon,**sigma, **lambda; + double nlambda, alphalj, alphac; + double **lj1,**lj2,**lj3,**lj4,**offset; + double *cut_respa; + double qdist; // TIP4P distance O to negative charge (compatibility of cutoffs) + double g_ewald; + + 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-FEP/pair_lj_cut_soft.cpp b/src/USER-FEP/pair_lj_cut_soft.cpp new file mode 100644 index 0000000000..22af520a2a --- /dev/null +++ b/src/USER-FEP/pair_lj_cut_soft.cpp @@ -0,0 +1,779 @@ +/* ---------------------------------------------------------------------- + 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) + Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lj_cut_soft.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "update.h" +#include "integrate.h" +#include "respa.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairLJCutSoft::PairLJCutSoft(LAMMPS *lmp) : Pair(lmp) +{ + respa_enable = 1; + writedata = 1; + allocated = 0; +} + +/* ---------------------------------------------------------------------- */ + +PairLJCutSoft::~PairLJCutSoft() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(epsilon); + memory->destroy(sigma); + memory->destroy(lambda); + memory->destroy(lj1); + memory->destroy(lj2); + memory->destroy(lj3); + memory->destroy(offset); + allocated=0; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutSoft::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,forcelj,factor_lj; + double denlj, r4sig6; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + 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]; + 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_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + + fpair = factor_lj*forcelj; + + 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) { + evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * + (1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype]; + evdwl *= factor_lj; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutSoft::compute_inner() +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fpair; + double rsq,forcelj,factor_lj,rsw; + double denlj, r4sig6; + int *ilist,*jlist,*numneigh,**firstneigh; + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + inum = listinner->inum; + ilist = listinner->ilist; + numneigh = listinner->numneigh; + firstneigh = listinner->firstneigh; + + double cut_out_on = cut_respa[0]; + double cut_out_off = cut_respa[1]; + + double cut_out_diff = cut_out_off - cut_out_on; + double cut_out_on_sq = cut_out_on*cut_out_on; + double cut_out_off_sq = cut_out_off*cut_out_off; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cut_out_off_sq) { + jtype = type[j]; + + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + + fpair = factor_lj*forcelj; + + if (rsq > cut_out_on_sq) { + rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; + fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw); + } + + 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; + } + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutSoft::compute_middle() +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fpair; + double rsq,forcelj,factor_lj,rsw; + double denlj, r4sig6; + int *ilist,*jlist,*numneigh,**firstneigh; + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + inum = listmiddle->inum; + ilist = listmiddle->ilist; + numneigh = listmiddle->numneigh; + firstneigh = listmiddle->firstneigh; + + double cut_in_off = cut_respa[0]; + double cut_in_on = cut_respa[1]; + double cut_out_on = cut_respa[2]; + double cut_out_off = cut_respa[3]; + + double cut_in_diff = cut_in_on - cut_in_off; + double cut_out_diff = cut_out_off - cut_out_on; + double cut_in_off_sq = cut_in_off*cut_in_off; + double cut_in_on_sq = cut_in_on*cut_in_on; + double cut_out_on_sq = cut_out_on*cut_out_on; + double cut_out_off_sq = cut_out_off*cut_out_off; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { + jtype = type[j]; + + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + + fpair = factor_lj*forcelj; + + if (rsq < cut_in_on_sq) { + rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; + fpair *= rsw*rsw*(3.0 - 2.0*rsw); + } + if (rsq > cut_out_on_sq) { + rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; + fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); + } + + 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; + } + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutSoft::compute_outer(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,forcelj,factor_lj,rsw; + double denlj, r4sig6; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + inum = listouter->inum; + ilist = listouter->ilist; + numneigh = listouter->numneigh; + firstneigh = listouter->firstneigh; + + double cut_in_off = cut_respa[2]; + double cut_in_on = cut_respa[3]; + + double cut_in_diff = cut_in_on - cut_in_off; + double cut_in_off_sq = cut_in_off*cut_in_off; + double cut_in_on_sq = cut_in_on*cut_in_on; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + if (rsq > cut_in_off_sq) { + + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + + fpair = factor_lj*forcelj; + + if (rsq < cut_in_on_sq) { + rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; + fpair *= rsw*rsw*(3.0 - 2.0*rsw); + } + + 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) { + evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * + (1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype]; + evdwl *= factor_lj; + } + + if (vflag) { + if (rsq <= cut_in_off_sq) { + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + fpair = factor_lj*forcelj; + } else if (rsq < cut_in_on_sq) + fpair = factor_lj*forcelj; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); + } + } + } +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairLJCutSoft::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(epsilon,n+1,n+1,"pair:epsilon"); + memory->create(sigma,n+1,n+1,"pair:sigma"); + memory->create(lambda,n+1,n+1,"pair:lambda"); + 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(offset,n+1,n+1,"pair:offset"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLJCutSoft::settings(int narg, char **arg) +{ + if (narg != 3) error->all(FLERR,"Illegal pair_style command"); + + nlambda = force->numeric(FLERR,arg[0]); + alphalj = force->numeric(FLERR,arg[1]); + + cut_global = 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]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJCutSoft::coeff(int narg, char **arg) +{ + if (narg < 5 || narg > 6) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double epsilon_one = force->numeric(FLERR,arg[2]); + double sigma_one = force->numeric(FLERR,arg[3]); + double lambda_one = force->numeric(FLERR,arg[4]); + + double cut_one = cut_global; + if (narg == 6) cut_one = force->numeric(FLERR,arg[5]); + + 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; + lambda[i][j] = lambda_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJCutSoft::init_style() +{ + // request regular or rRESPA neighbor lists + + int irequest; + + if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { + int respa = 0; + if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; + if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + + if (respa == 0) irequest = neighbor->request(this); + else if (respa == 1) { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } else { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 2; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respamiddle = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } + + } else irequest = neighbor->request(this); + + // set rRESPA cutoffs + + if (strstr(update->integrate_style,"respa") && + ((Respa *) update->integrate)->level_inner >= 0) + cut_respa = ((Respa *) update->integrate)->cutoff; + else cut_respa = NULL; + +} + +/* ---------------------------------------------------------------------- + neighbor callback to inform pair style of neighbor list to use + regular or rRESPA +------------------------------------------------------------------------- */ + +void PairLJCutSoft::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 PairLJCutSoft::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[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + lj1[i][j] = pow(lambda[i][j], nlambda); + lj2[i][j] = pow(sigma[i][j], 6.0); + lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]); + + if (offset_flag) { + double denlj = lj3[i][j] + pow(sigma[i][j] / cut[i][j], 6.0); + offset[i][j] = lj1[i][j] * 4.0 * epsilon[i][j] * (1.0/(denlj*denlj) - 1.0/denlj); + } else offset[i][j] = 0.0; + + epsilon[j][i] = epsilon[i][j]; + sigma[j][i] = sigma[i][j]; + lambda[j][i] = lambda[i][j]; + lj1[j][i] = lj1[i][j]; + lj2[j][i] = lj2[i][j]; + lj3[j][i] = lj3[i][j]; + offset[j][i] = offset[i][j]; + + // check interior rRESPA cutoff + + if (cut_respa && cut[i][j] < 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[i][j]*cut[i][j]*cut[i][j]; + double rc6 = rc3*rc3; + double rc9 = rc3*rc6; + etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * + sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); + ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * + sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); + } + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJCutSoft::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(&lambda[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutSoft::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(&lambda[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&lambda[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJCutSoft::write_restart_settings(FILE *fp) +{ + fwrite(&nlambda,sizeof(double),1,fp); + fwrite(&alphalj,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); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutSoft::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + fread(&nlambda,sizeof(double),1,fp); + fread(&alphalj,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); + } + MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world); + MPI_Bcast(&alphalj,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); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairLJCutSoft::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],lambda[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairLJCutSoft::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 %g\n",i,j,epsilon[i][j],sigma[i][j], + lambda[i][j],cut[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJCutSoft::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double forcelj,philj; + double r4sig6, denlj; + + if (rsq < cutsq[itype][jtype]) { + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + } else forcelj = 0.0; + fforce = factor_lj*forcelj; + + if (rsq < cutsq[itype][jtype]) { + philj = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * + (1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype]; + } else philj = 0.0; + + return factor_lj*philj; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJCutSoft::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str,"lambda") == 0) return (void *) lambda; + return NULL; +} diff --git a/src/USER-FEP/pair_lj_cut_soft.h b/src/USER-FEP/pair_lj_cut_soft.h new file mode 100644 index 0000000000..50ce685e5c --- /dev/null +++ b/src/USER-FEP/pair_lj_cut_soft.h @@ -0,0 +1,87 @@ +/* -*- 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/soft,PairLJCutSoft) + +#else + +#ifndef LMP_PAIR_LJ_CUT_SOFT_H +#define LMP_PAIR_LJ_CUT_SOFT_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLJCutSoft : public Pair { + public: + PairLJCutSoft(class LAMMPS *); + virtual ~PairLJCutSoft(); + 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 &); + void *extract(const char *, int &); + + void compute_inner(); + void compute_middle(); + void compute_outer(int, int); + + protected: + double cut_global; + double **cut; + double **epsilon,**sigma, **lambda; + double nlambda, alphalj; + double **lj1,**lj2,**lj3,**offset; + double *cut_respa; + + 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 cutoff < Respa interior cutoff + +One or more pairwise cutoffs are too short to use with the specified +rRESPA cutoffs. + +E: Pair lj/cut/soft different lambda values in mix + +The value of lambda has to be the same for I J pairs. + +*/ diff --git a/src/USER-FEP/pair_lj_cut_tip4p_long_soft.cpp b/src/USER-FEP/pair_lj_cut_tip4p_long_soft.cpp new file mode 100644 index 0000000000..0ea6cf7a7f --- /dev/null +++ b/src/USER-FEP/pair_lj_cut_tip4p_long_soft.cpp @@ -0,0 +1,592 @@ +/* ---------------------------------------------------------------------- + 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 authors: Amalie Frischknecht and Ahmed Ismail (SNL) + simpler force assignment added by Rolf Isele-Holder (Aachen University) + Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lj_cut_tip4p_long_soft.h" +#include "angle.h" +#include "atom.h" +#include "bond.h" +#include "comm.h" +#include "domain.h" +#include "force.h" +#include "kspace.h" +#include "update.h" +#include "respa.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairLJCutTIP4PLongSoft::PairLJCutTIP4PLongSoft(LAMMPS *lmp) : + PairLJCutCoulLongSoft(lmp) +{ + tip4pflag = 1; + single_enable = 0; + respa_enable = 0; + + nmax = 0; + hneigh = NULL; + newsite = NULL; + + // TIP4P cannot compute virial as F dot r + // due to finding bonded H atoms which are not near O atom + + no_virial_fdotr_compute = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairLJCutTIP4PLongSoft::~PairLJCutTIP4PLongSoft() +{ + memory->destroy(hneigh); + memory->destroy(newsite); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype,key; + int n,vlist[6]; + int iH1,iH2,jH1,jH2; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul; + double r,forcecoul,forcelj,cforce; + double factor_coul,factor_lj; + double grij,expm2,prefactor,t,erfc; + double denc, denlj, r4sig6; + double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3]; + double *x1,*x2; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + // reallocate hneigh & newsite if necessary + // initialize hneigh[0] to -1 on steps when reneighboring occurred + // initialize hneigh[2] to 0 every step + + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + if (atom->nmax > nmax) { + nmax = atom->nmax; + memory->destroy(hneigh); + memory->create(hneigh,nmax,3,"pair:hneigh"); + memory->destroy(newsite); + memory->create(newsite,nmax,3,"pair:newsite"); + } + if (neighbor->ago == 0) + for (i = 0; i < nall; i++) hneigh[i][0] = -1; + for (i = 0; i < nall; i++) hneigh[i][2] = 0; + + double **f = atom->f; + double **x = atom->x; + double *q = atom->q; + int *type = atom->type; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + double cut_coulsqplus = (cut_coul+2.0*qdist) * (cut_coul+2.0*qdist); + + 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]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + + // if atom I = water O, set x1 = offset charge site + // else x1 = x of atom I + + if (itype == typeO) { + if (hneigh[i][0] < 0) { + hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1); + hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2); + hneigh[i][2] = 1; + if (iH1 == -1 || iH2 == -1) + error->one(FLERR,"TIP4P hydrogen is missing"); + if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) + error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); + compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); + } else { + iH1 = hneigh[i][0]; + iH2 = hneigh[i][1]; + if (hneigh[i][2] == 0) { + hneigh[i][2] = 1; + compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); + } + } + x1 = newsite[i]; + } else x1 = x[i]; + + jlist = firstneigh[i]; + jnum = numneigh[i]; + + 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]; + + // LJ interaction based on true rsq + + if (rsq < cut_ljsq[itype][jtype]) { + + r4sig6 = rsq*rsq / lj2[itype][jtype]; + denlj = lj3[itype][jtype] + rsq*r4sig6; + forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * + (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj)); + + forcelj *= factor_lj; + + f[i][0] += delx*forcelj; + f[i][1] += dely*forcelj; + f[i][2] += delz*forcelj; + f[j][0] -= delx*forcelj; + f[j][1] -= dely*forcelj; + f[j][2] -= delz*forcelj; + + if (eflag) { + evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * + (1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,forcelj,delx,dely,delz); + } + + // adjust rsq and delxyz for off-site O charge(s) if necessary + // but only if they are within reach + + if (rsq < cut_coulsqplus) { + if (itype == typeO || jtype == typeO) { + + // if atom J = water O, set x2 = offset charge site + // else x2 = x of atom J + + if (jtype == typeO) { + if (hneigh[j][0] < 0) { + hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1); + hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2); + hneigh[j][2] = 1; + if (jH1 == -1 || jH2 == -1) + error->one(FLERR,"TIP4P hydrogen is missing"); + if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) + error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); + compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); + } else { + jH1 = hneigh[j][0]; + jH2 = hneigh[j][1]; + if (hneigh[j][2] == 0) { + hneigh[j][2] = 1; + compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); + } + } + x2 = newsite[j]; + } else x2 = x[j]; + + delx = x1[0] - x2[0]; + dely = x1[1] - x2[1]; + delz = x1[2] - x2[2]; + rsq = delx*delx + dely*dely + delz*delz; + } + + // Coulombic interaction based on modified rsq + + if (rsq < cut_coulsq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + + denc = sqrt(lj4[itype][jtype] + rsq); + prefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc); + + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) { + forcecoul -= (1.0-factor_coul)*prefactor; + } + + cforce = forcecoul; + + // if i,j are not O atoms, force is applied directly + // if i or j are O atoms, force is on fictitious atom & partitioned + // force partitioning due to Feenstra, J Comp Chem, 20, 786 (1999) + // f_f = fictitious force, fO = f_f (1 - 2 alpha), fH = alpha f_f + // preserves total force and torque on water molecule + // virial = sum(r x F) where each water's atoms are near xi and xj + // vlist stores 2,4,6 atoms whose forces contribute to virial + + n = 0; + key = 0; + + if (itype != typeO) { + f[i][0] += delx * cforce; + f[i][1] += dely * cforce; + f[i][2] += delz * cforce; + + if (vflag) { + v[0] = x[i][0] * delx * cforce; + v[1] = x[i][1] * dely * cforce; + v[2] = x[i][2] * delz * cforce; + v[3] = x[i][0] * dely * cforce; + v[4] = x[i][0] * delz * cforce; + v[5] = x[i][1] * delz * cforce; + } + vlist[n++] = i; + + } else { + key++; + + fd[0] = delx*cforce; + fd[1] = dely*cforce; + fd[2] = delz*cforce; + + fO[0] = fd[0]*(1 - alpha); + fO[1] = fd[1]*(1 - alpha); + fO[2] = fd[2]*(1 - alpha); + + fH[0] = 0.5 * alpha * fd[0]; + fH[1] = 0.5 * alpha * fd[1]; + fH[2] = 0.5 * alpha * fd[2]; + + f[i][0] += fO[0]; + f[i][1] += fO[1]; + f[i][2] += fO[2]; + + f[iH1][0] += fH[0]; + f[iH1][1] += fH[1]; + f[iH1][2] += fH[2]; + + f[iH2][0] += fH[0]; + f[iH2][1] += fH[1]; + f[iH2][2] += fH[2]; + + if (vflag) { + domain->closest_image(x[i],x[iH1],xH1); + domain->closest_image(x[i],x[iH2],xH2); + + v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; + v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; + v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; + v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; + v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; + v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2]; + } + vlist[n++] = i; + vlist[n++] = iH1; + vlist[n++] = iH2; + } + + if (jtype != typeO) { + f[j][0] -= delx * cforce; + f[j][1] -= dely * cforce; + f[j][2] -= delz * cforce; + + if (vflag) { + v[0] -= x[j][0] * delx * cforce; + v[1] -= x[j][1] * dely * cforce; + v[2] -= x[j][2] * delz * cforce; + v[3] -= x[j][0] * dely * cforce; + v[4] -= x[j][0] * delz * cforce; + v[5] -= x[j][1] * delz * cforce; + } + vlist[n++] = j; + + } else { + key += 2; + + fd[0] = -delx*cforce; + fd[1] = -dely*cforce; + fd[2] = -delz*cforce; + + fO[0] = fd[0]*(1 - alpha); + fO[1] = fd[1]*(1 - alpha); + fO[2] = fd[2]*(1 - alpha); + + fH[0] = 0.5 * alpha * fd[0]; + fH[1] = 0.5 * alpha * fd[1]; + fH[2] = 0.5 * alpha * fd[2]; + + f[j][0] += fO[0]; + f[j][1] += fO[1]; + f[j][2] += fO[2]; + + f[jH1][0] += fH[0]; + f[jH1][1] += fH[1]; + f[jH1][2] += fH[2]; + + f[jH2][0] += fH[0]; + f[jH2][1] += fH[1]; + f[jH2][2] += fH[2]; + + if (vflag) { + domain->closest_image(x[j],x[jH1],xH1); + domain->closest_image(x[j],x[jH2],xH2); + + v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; + v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; + v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; + v[3] += x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; + v[4] += x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; + v[5] += x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2]; + } + vlist[n++] = j; + vlist[n++] = jH1; + vlist[n++] = jH2; + } + + if (eflag) { + prefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / denc; + ecoul = prefactor*erfc; + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + + if (evflag) ev_tally_tip4p(key,vlist,v,ecoul,alpha); + } + } + } + } +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLJCutTIP4PLongSoft::settings(int narg, char **arg) +{ + if (narg < 9 || narg > 10) error->all(FLERR,"Illegal pair_style command"); + + typeO = force->inumeric(FLERR,arg[0]); + typeH = force->inumeric(FLERR,arg[1]); + typeB = force->inumeric(FLERR,arg[2]); + typeA = force->inumeric(FLERR,arg[3]); + qdist = force->numeric(FLERR,arg[4]); + nlambda = force->numeric(FLERR,arg[5]); + alphalj = force->numeric(FLERR,arg[6]); + alphac = force->numeric(FLERR,arg[7]); + + cut_lj_global = force->numeric(FLERR,arg[8]); + if (narg == 9) cut_coul = cut_lj_global; + else cut_coul = force->numeric(FLERR,arg[9]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; + } +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJCutTIP4PLongSoft::init_style() +{ + if (atom->tag_enable == 0) + error->all(FLERR,"Pair style lj/cut/tip4p/long/soft requires atom IDs"); + if (!force->newton_pair) + error->all(FLERR, + "Pair style lj/cut/tip4p/long/soft requires newton pair on"); + if (!atom->q_flag) + error->all(FLERR, + "Pair style lj/cut/tip4p/long/soft requires atom attribute q"); + if (force->bond == NULL) + error->all(FLERR,"Must use a bond style with TIP4P potential"); + if (force->angle == NULL) + error->all(FLERR,"Must use an angle style with TIP4P potential"); + + PairLJCutCoulLongSoft::init_style(); + + // set alpha parameter + + double theta = force->angle->equilibrium_angle(typeA); + double blen = force->bond->equilibrium_distance(typeB); + alpha = qdist / (cos(0.5*theta) * blen); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLJCutTIP4PLongSoft::init_one(int i, int j) +{ + double cut = PairLJCutCoulLongSoft::init_one(i,j); + + // check that LJ epsilon = 0.0 for water H + // set LJ cutoff to 0.0 for any interaction involving water H + // so LJ term isn't calculated in compute() + + if ((i == typeH && epsilon[i][i] != 0.0) || + (j == typeH && epsilon[j][j] != 0.0)) + error->all(FLERR,"Water H epsilon must be 0.0 for " + "pair style lj/cut/tip4p/long/soft"); + + if (i == typeH || j == typeH) + cut_ljsq[j][i] = cut_ljsq[i][j] = 0.0; + + return cut; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJCutTIP4PLongSoft::write_restart_settings(FILE *fp) +{ + PairLJCutCoulLongSoft::write_restart_settings(fp); + + fwrite(&typeO,sizeof(int),1,fp); + fwrite(&typeH,sizeof(int),1,fp); + fwrite(&typeB,sizeof(int),1,fp); + fwrite(&typeA,sizeof(int),1,fp); + fwrite(&qdist,sizeof(double),1,fp); + + fwrite(&cut_lj_global,sizeof(double),1,fp); + fwrite(&cut_coul,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); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutTIP4PLongSoft::read_restart_settings(FILE *fp) +{ + PairLJCutCoulLongSoft::read_restart_settings(fp); + + if (comm->me == 0) { + fread(&typeO,sizeof(int),1,fp); + fread(&typeH,sizeof(int),1,fp); + fread(&typeB,sizeof(int),1,fp); + fread(&typeA,sizeof(int),1,fp); + fread(&qdist,sizeof(double),1,fp); + + fread(&cut_lj_global,sizeof(double),1,fp); + fread(&cut_coul,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); + } + + MPI_Bcast(&typeO,1,MPI_INT,0,world); + MPI_Bcast(&typeH,1,MPI_INT,0,world); + MPI_Bcast(&typeB,1,MPI_INT,0,world); + MPI_Bcast(&typeA,1,MPI_INT,0,world); + MPI_Bcast(&qdist,1,MPI_DOUBLE,0,world); + + MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul,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); +} + +/* ---------------------------------------------------------------------- + compute position xM of fictitious charge site for O atom and 2 H atoms + return it as xM +------------------------------------------------------------------------- */ + +void PairLJCutTIP4PLongSoft::compute_newsite(double *xO, double *xH1, + double *xH2, double *xM) +{ + double delx1 = xH1[0] - xO[0]; + double dely1 = xH1[1] - xO[1]; + double delz1 = xH1[2] - xO[2]; + domain->minimum_image(delx1,dely1,delz1); + + double delx2 = xH2[0] - xO[0]; + double dely2 = xH2[1] - xO[1]; + double delz2 = xH2[2] - xO[2]; + domain->minimum_image(delx2,dely2,delz2); + + xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2); + xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2); + xM[2] = xO[2] + alpha * 0.5 * (delz1 + delz2); +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJCutTIP4PLongSoft::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"qdist") == 0) return (void *) &qdist; + if (strcmp(str,"typeO") == 0) return (void *) &typeO; + if (strcmp(str,"typeH") == 0) return (void *) &typeH; + if (strcmp(str,"typeA") == 0) return (void *) &typeA; + if (strcmp(str,"typeB") == 0) return (void *) &typeB; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str,"lambda") == 0) return (void *) lambda; + return NULL; +} + +/* ---------------------------------------------------------------------- + memory usage of hneigh +------------------------------------------------------------------------- */ + +double PairLJCutTIP4PLongSoft::memory_usage() +{ + double bytes = maxeatom * sizeof(double); + bytes += maxvatom*6 * sizeof(double); + bytes += 2 * nmax * sizeof(double); + return bytes; +} diff --git a/src/USER-FEP/pair_lj_cut_tip4p_long_soft.h b/src/USER-FEP/pair_lj_cut_tip4p_long_soft.h new file mode 100644 index 0000000000..5255882a00 --- /dev/null +++ b/src/USER-FEP/pair_lj_cut_tip4p_long_soft.h @@ -0,0 +1,105 @@ +/* ---------------------------------------------------------------------- + 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/tip4p/long/soft,PairLJCutTIP4PLongSoft) + +#else + +#ifndef LMP_PAIR_LJ_CUT_TIP4P_LONG_SOFT_H +#define LMP_PAIR_LJ_CUT_TIP4P_LONG_SOFT_H + +#include "pair_lj_cut_coul_long_soft.h" + +namespace LAMMPS_NS { + +class PairLJCutTIP4PLongSoft : public PairLJCutCoulLongSoft { + public: + PairLJCutTIP4PLongSoft(class LAMMPS *); + virtual ~PairLJCutTIP4PLongSoft(); + virtual void compute(int, int); + void settings(int, char **); + void init_style(); + double init_one(int, int); + void write_restart_settings(FILE *fp); + void read_restart_settings(FILE *fp); + void *extract(const char *, int &); + virtual double memory_usage(); + + protected: + int typeH,typeO; // atom types of TIP4P water H and O atoms + int typeA,typeB; // angle and bond types of TIP4P water + double alpha; // geometric constraint parameter for TIP4P + + int nmax; // info on off-oxygen charge sites + int **hneigh; // 0,1 = indices of 2 H associated with O + // 2 = 0 if site loc not yet computed, 1 if yes + double **newsite; // locations of charge sites + + void compute_newsite(double *, double *, double *, double *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: TIP4P hydrogen is missing + +The TIP4P pairwise computation failed to find the correct H atom +within a water molecule. + +E: TIP4P hydrogen has incorrect atom type + +The TIP4P pairwise computation found an H atom whose type does not +agree with the specified H type. + +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: Pair style lj/cut/tip4p/long/soft requires atom IDs + +There are no atom IDs defined in the system and the TIP4P potential +requires them to find O,H atoms with a water molecule. + +E: Pair style lj/cut/tip4p/long/soft requires newton pair on + +This is because the computation of constraint forces within a water +molecule adds forces to atoms owned by other processors. + +E: Pair style lj/cut/tip4p/long/soft requires atom attribute q + +The atom style defined does not have these attributes. + +E: Must use a bond style with TIP4P potential + +TIP4P potentials assume bond lengths in water are constrained +by a fix shake command. + +E: Must use an angle style with TIP4P potential + +TIP4P potentials assume angles in water are constrained by a fix shake +command. + +E: Water H epsilon must be 0.0 for pair style lj/cut/tip4p/long/soft + +This is because LAMMPS does not compute the Lennard-Jones interactions +with these particles for efficiency reasons. + +*/ diff --git a/src/USER-FEP/pair_tip4p_long_soft.cpp b/src/USER-FEP/pair_tip4p_long_soft.cpp new file mode 100644 index 0000000000..ef76da4994 --- /dev/null +++ b/src/USER-FEP/pair_tip4p_long_soft.cpp @@ -0,0 +1,514 @@ +/* ---------------------------------------------------------------------- + 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 authors: Amalie Frischknecht and Ahmed Ismail (SNL) + simpler force assignment added by Rolf Isele-Holder (Aachen University) + Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_tip4p_long_soft.h" +#include "angle.h" +#include "atom.h" +#include "bond.h" +#include "comm.h" +#include "domain.h" +#include "force.h" +#include "kspace.h" +#include "update.h" +#include "respa.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairTIP4PLongSoft::PairTIP4PLongSoft(LAMMPS *lmp) : PairCoulLongSoft(lmp) +{ + tip4pflag = 1; + single_enable = 0; + respa_enable = 0; + + nmax = 0; + hneigh = NULL; + newsite = NULL; + + // TIP4P cannot compute virial as F dot r + // due to finding bonded H atoms which are not near O atom + + no_virial_fdotr_compute = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairTIP4PLongSoft::~PairTIP4PLongSoft() +{ + memory->destroy(hneigh); + memory->destroy(newsite); +} + +/* ---------------------------------------------------------------------- */ + +void PairTIP4PLongSoft::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype,key; + int n,vlist[6]; + int iH1,iH2,jH1,jH2; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul; + double r,forcecoul,cforce; + double factor_coul; + double grij,expm2,prefactor,t,erfc; + double denc; + double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3]; + double *x1,*x2; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + // reallocate hneigh & newsite if necessary + // initialize hneigh[0] to -1 on steps when reneighboring occurred + // initialize hneigh[2] to 0 every step + + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + if (atom->nmax > nmax) { + nmax = atom->nmax; + memory->destroy(hneigh); + memory->create(hneigh,nmax,3,"pair:hneigh"); + memory->destroy(newsite); + memory->create(newsite,nmax,3,"pair:newsite"); + } + if (neighbor->ago == 0) + for (i = 0; i < nall; i++) hneigh[i][0] = -1; + for (i = 0; i < nall; i++) hneigh[i][2] = 0; + + double **f = atom->f; + double **x = atom->x; + double *q = atom->q; + int *type = atom->type; + double *special_coul = force->special_coul; + double qqrd2e = force->qqrd2e; + double cut_coulsqplus = (cut_coul+2.0*qdist) * (cut_coul+2.0*qdist); + + 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]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + + // if atom I = water O, set x1 = offset charge site + // else x1 = x of atom I + + if (itype == typeO) { + if (hneigh[i][0] < 0) { + hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1); + hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2); + hneigh[i][2] = 1; + if (iH1 == -1 || iH2 == -1) + error->one(FLERR,"TIP4P hydrogen is missing"); + if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) + error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); + compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); + } else { + iH1 = hneigh[i][0]; + iH2 = hneigh[i][1]; + if (hneigh[i][2] == 0) { + hneigh[i][2] = 1; + compute_newsite(x[i],x[iH1],x[iH2],newsite[i]); + } + } + x1 = newsite[i]; + } else x1 = x[i]; + + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + 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]; + + // adjust rsq and delxyz for off-site O charge(s) if necessary + // but only if they are within reach + + if (rsq < cut_coulsqplus) { + if (itype == typeO || jtype == typeO) { + + // if atom J = water O, set x2 = offset charge site + // else x2 = x of atom J + + if (jtype == typeO) { + if (hneigh[j][0] < 0) { + hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1); + hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2); + hneigh[j][2] = 1; + if (jH1 == -1 || jH2 == -1) + error->one(FLERR,"TIP4P hydrogen is missing"); + if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) + error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); + compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); + } else { + jH1 = hneigh[j][0]; + jH2 = hneigh[j][1]; + if (hneigh[j][2] == 0) { + hneigh[j][2] = 1; + compute_newsite(x[j],x[jH1],x[jH2],newsite[j]); + } + } + x2 = newsite[j]; + } else x2 = x[j]; + + delx = x1[0] - x2[0]; + dely = x1[1] - x2[1]; + delz = x1[2] - x2[2]; + rsq = delx*delx + dely*dely + delz*delz; + } + + // Coulombic interaction based on modified rsq + + if (rsq < cut_coulsq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + + denc = sqrt(lam2[itype][jtype] + rsq); + prefactor = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / (denc*denc*denc); + + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) { + forcecoul -= (1.0-factor_coul)*prefactor; + } + + cforce = forcecoul; + + // if i,j are not O atoms, force is applied directly + // if i or j are O atoms, force is on fictitious atom & partitioned + // force partitioning due to Feenstra, J Comp Chem, 20, 786 (1999) + // f_f = fictitious force, fO = f_f (1 - 2 alpha), fH = alpha f_f + // preserves total force and torque on water molecule + // virial = sum(r x F) where each water's atoms are near xi and xj + // vlist stores 2,4,6 atoms whose forces contribute to virial + + n = 0; + key = 0; + + if (itype != typeO) { + f[i][0] += delx * cforce; + f[i][1] += dely * cforce; + f[i][2] += delz * cforce; + + if (vflag) { + v[0] = x[i][0] * delx * cforce; + v[1] = x[i][1] * dely * cforce; + v[2] = x[i][2] * delz * cforce; + v[3] = x[i][0] * dely * cforce; + v[4] = x[i][0] * delz * cforce; + v[5] = x[i][1] * delz * cforce; + } + vlist[n++] = i; + + } else { + key++; + + fd[0] = delx*cforce; + fd[1] = dely*cforce; + fd[2] = delz*cforce; + + fO[0] = fd[0]*(1 - alpha); + fO[1] = fd[1]*(1 - alpha); + fO[2] = fd[2]*(1 - alpha); + + fH[0] = 0.5 * alpha * fd[0]; + fH[1] = 0.5 * alpha * fd[1]; + fH[2] = 0.5 * alpha * fd[2]; + + f[i][0] += fO[0]; + f[i][1] += fO[1]; + f[i][2] += fO[2]; + + f[iH1][0] += fH[0]; + f[iH1][1] += fH[1]; + f[iH1][2] += fH[2]; + + f[iH2][0] += fH[0]; + f[iH2][1] += fH[1]; + f[iH2][2] += fH[2]; + + if (vflag) { + domain->closest_image(x[i],x[iH1],xH1); + domain->closest_image(x[i],x[iH2],xH2); + + v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; + v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; + v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; + v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; + v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; + v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2]; + } + vlist[n++] = i; + vlist[n++] = iH1; + vlist[n++] = iH2; + } + + if (jtype != typeO) { + f[j][0] -= delx * cforce; + f[j][1] -= dely * cforce; + f[j][2] -= delz * cforce; + + if (vflag) { + v[0] -= x[j][0] * delx * cforce; + v[1] -= x[j][1] * dely * cforce; + v[2] -= x[j][2] * delz * cforce; + v[3] -= x[j][0] * dely * cforce; + v[4] -= x[j][0] * delz * cforce; + v[5] -= x[j][1] * delz * cforce; + } + vlist[n++] = j; + + } else { + key += 2; + + fd[0] = -delx*cforce; + fd[1] = -dely*cforce; + fd[2] = -delz*cforce; + + fO[0] = fd[0]*(1 - alpha); + fO[1] = fd[1]*(1 - alpha); + fO[2] = fd[2]*(1 - alpha); + + fH[0] = 0.5 * alpha * fd[0]; + fH[1] = 0.5 * alpha * fd[1]; + fH[2] = 0.5 * alpha * fd[2]; + + f[j][0] += fO[0]; + f[j][1] += fO[1]; + f[j][2] += fO[2]; + + f[jH1][0] += fH[0]; + f[jH1][1] += fH[1]; + f[jH1][2] += fH[2]; + + f[jH2][0] += fH[0]; + f[jH2][1] += fH[1]; + f[jH2][2] += fH[2]; + + if (vflag) { + domain->closest_image(x[j],x[jH1],xH1); + domain->closest_image(x[j],x[jH2],xH2); + + v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; + v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; + v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; + v[3] += x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; + v[4] += x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; + v[5] += x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2]; + } + vlist[n++] = j; + vlist[n++] = jH1; + vlist[n++] = jH2; + } + + if (eflag) { + prefactor = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / denc; + ecoul = prefactor*erfc; + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + + if (evflag) ev_tally_tip4p(key,vlist,v,ecoul,alpha); + } + } + } + } +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairTIP4PLongSoft::settings(int narg, char **arg) +{ + if (narg != 8) error->all(FLERR,"Illegal pair_style command"); + + typeO = force->inumeric(FLERR,arg[0]); + typeH = force->inumeric(FLERR,arg[1]); + typeB = force->inumeric(FLERR,arg[2]); + typeA = force->inumeric(FLERR,arg[3]); + qdist = force->numeric(FLERR,arg[4]); + + nlambda = force->numeric(FLERR,arg[5]); + alphac = force->numeric(FLERR,arg[6]); + + cut_coul = force->numeric(FLERR,arg[7]); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairTIP4PLongSoft::init_style() +{ + if (atom->tag_enable == 0) + error->all(FLERR,"Pair style tip4p/long requires atom IDs"); + if (!force->newton_pair) + error->all(FLERR, + "Pair style tip4p/long requires newton pair on"); + if (!atom->q_flag) + error->all(FLERR, + "Pair style tip4p/long requires atom attribute q"); + if (force->bond == NULL) + error->all(FLERR,"Must use a bond style with TIP4P potential"); + if (force->angle == NULL) + error->all(FLERR,"Must use an angle style with TIP4P potential"); + + PairCoulLongSoft::init_style(); + + // set alpha parameter + + double theta = force->angle->equilibrium_angle(typeA); + double blen = force->bond->equilibrium_distance(typeB); + alpha = qdist / (cos(0.5*theta) * blen); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairTIP4PLongSoft::init_one(int i, int j) +{ + return PairCoulLongSoft::init_one(i,j); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairTIP4PLongSoft::write_restart_settings(FILE *fp) +{ + PairCoulLongSoft::write_restart_settings(fp); + + fwrite(&typeO,sizeof(int),1,fp); + fwrite(&typeH,sizeof(int),1,fp); + fwrite(&typeB,sizeof(int),1,fp); + fwrite(&typeA,sizeof(int),1,fp); + fwrite(&qdist,sizeof(double),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairTIP4PLongSoft::read_restart_settings(FILE *fp) +{ + PairCoulLongSoft::read_restart_settings(fp); + + if (comm->me == 0) { + fread(&typeO,sizeof(int),1,fp); + fread(&typeH,sizeof(int),1,fp); + fread(&typeB,sizeof(int),1,fp); + fread(&typeA,sizeof(int),1,fp); + fread(&qdist,sizeof(double),1,fp); + } + MPI_Bcast(&typeO,1,MPI_INT,0,world); + MPI_Bcast(&typeH,1,MPI_INT,0,world); + MPI_Bcast(&typeB,1,MPI_INT,0,world); + MPI_Bcast(&typeA,1,MPI_INT,0,world); + MPI_Bcast(&qdist,1,MPI_DOUBLE,0,world); +} + +/* ---------------------------------------------------------------------- + compute position xM of fictitious charge site for O atom and 2 H atoms + return it as xM +------------------------------------------------------------------------- */ + +void PairTIP4PLongSoft::compute_newsite(double *xO, double *xH1, + double *xH2, double *xM) +{ + double delx1 = xH1[0] - xO[0]; + double dely1 = xH1[1] - xO[1]; + double delz1 = xH1[2] - xO[2]; + domain->minimum_image(delx1,dely1,delz1); + + double delx2 = xH2[0] - xO[0]; + double dely2 = xH2[1] - xO[1]; + double delz2 = xH2[2] - xO[2]; + domain->minimum_image(delx2,dely2,delz2); + + xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2); + xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2); + xM[2] = xO[2] + alpha * 0.5 * (delz1 + delz2); +} + +/* ---------------------------------------------------------------------- */ + +void *PairTIP4PLongSoft::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"qdist") == 0) return (void *) &qdist; + if (strcmp(str,"typeO") == 0) return (void *) &typeO; + if (strcmp(str,"typeH") == 0) return (void *) &typeH; + if (strcmp(str,"typeA") == 0) return (void *) &typeA; + if (strcmp(str,"typeB") == 0) return (void *) &typeB; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + dim = 2; + if (strcmp(str,"lambda") == 0) return (void *) lambda; + return NULL; +} + +/* ---------------------------------------------------------------------- + memory usage of hneigh +------------------------------------------------------------------------- */ + +double PairTIP4PLongSoft::memory_usage() +{ + double bytes = maxeatom * sizeof(double); + bytes += maxvatom*6 * sizeof(double); + bytes += 2 * nmax * sizeof(double); + return bytes; +} diff --git a/src/USER-FEP/pair_tip4p_long_soft.h b/src/USER-FEP/pair_tip4p_long_soft.h new file mode 100644 index 0000000000..c397cdbe21 --- /dev/null +++ b/src/USER-FEP/pair_tip4p_long_soft.h @@ -0,0 +1,100 @@ +/* ---------------------------------------------------------------------- + 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(tip4p/long/soft,PairTIP4PLongSoft) + +#else + +#ifndef LMP_PAIR_TIP4P_LONG_SOFT_H +#define LMP_PAIR_TIP4P_LONG_SOFT_H + +#include "pair_coul_long_soft.h" + +namespace LAMMPS_NS { + +class PairTIP4PLongSoft : public PairCoulLongSoft { + public: + PairTIP4PLongSoft(class LAMMPS *); + virtual ~PairTIP4PLongSoft(); + virtual void compute(int, int); + void settings(int, char **); + void init_style(); + double init_one(int, int); + void write_restart_settings(FILE *fp); + void read_restart_settings(FILE *fp); + void *extract(const char *, int &); + virtual double memory_usage(); + + protected: + int typeH,typeO; // atom types of TIP4P water H and O atoms + int typeA,typeB; // angle and bond types of TIP4P water + double alpha; // geometric constraint parameter for TIP4P + + int nmax; // info on off-oxygen charge sites + int **hneigh; // 0,1 = indices of 2 H associated with O + // 2 = 0 if site loc not yet computed, 1 if yes + double **newsite; // locations of charge sites + + void compute_newsite(double *, double *, double *, double *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: TIP4P hydrogen is missing + +The TIP4P pairwise computation failed to find the correct H atom +within a water molecule. + +E: TIP4P hydrogen has incorrect atom type + +The TIP4P pairwise computation found an H atom whose type does not +agree with the specified H type. + +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: Pair style tip4p/long/soft requires atom IDs + +There are no atom IDs defined in the system and the TIP4P potential +requires them to find O,H atoms with a water molecule. + +E: Pair style tip4p/long/soft requires newton pair on + +This is because the computation of constraint forces within a water +molecule adds forces to atoms owned by other processors. + +E: Pair style tip4p/long/soft requires atom attribute q + +The atom style defined does not have these attributes. + +E: Must use a bond style with TIP4P potential + +TIP4P potentials assume bond lengths in water are constrained +by a fix shake command. + +E: Must use an angle style with TIP4P potential + +TIP4P potentials assume angles in water are constrained by a fix shake +command. + +*/