diff --git a/src/USER-MISC/pair_lj_sf.cpp b/src/USER-MISC/pair_lj_sf.cpp index d378cfab4..2e1c92c93 100644 --- a/src/USER-MISC/pair_lj_sf.cpp +++ b/src/USER-MISC/pair_lj_sf.cpp @@ -1,351 +1,356 @@ /* ---------------------------------------------------------------------- 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: Laurent Joly (U Lyon), ljoly.ulyon@gmail.com ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "stdlib.h" #include "pair_lj_sf.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neigh_list.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairLJShiftedForce::PairLJShiftedForce(LAMMPS *lmp) : Pair(lmp) {} +PairLJShiftedForce::PairLJShiftedForce(LAMMPS *lmp) : Pair(lmp) +{ + respa_enable = 0; +} /* ---------------------------------------------------------------------- */ PairLJShiftedForce::~PairLJShiftedForce() { if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); memory->destroy(cut); memory->destroy(epsilon); memory->destroy(sigma); memory->destroy(lj1); memory->destroy(lj2); memory->destroy(lj3); memory->destroy(lj4); memory->destroy(foffset); memory->destroy(offset); } } /* ---------------------------------------------------------------------- */ void PairLJShiftedForce::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,r2inv,r6inv,forcelj,factor_lj; double r,t; 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; int nall = nlocal + atom->nghost; 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]) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; r = sqrt(rsq); t = r/cut[itype][jtype]; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) - t*foffset[itype][jtype]; fpair = factor_lj*forcelj*r2inv; f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } if (eflag) { evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) + (t-1.0)*foffset[itype][jtype] - 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(); } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairLJShiftedForce::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(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj3,n+1,n+1,"pair:lj2"); + 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(foffset,n+1,n+1,"pair:foffset"); memory->create(offset,n+1,n+1,"pair:offset"); } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairLJShiftedForce::settings(int narg, char **arg) { if (narg != 1) error->all(FLERR,"Illegal pair_style command"); cut_global = force->numeric(arg[0]); if (cut_global <= 0.0) error->all(FLERR,"Illegal pair_style command"); // 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 PairLJShiftedForce::coeff(int narg, char **arg) { - if (narg < 4 || narg > 5) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 4 || narg > 5) + error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); double epsilon_one = force->numeric(arg[2]); double sigma_one = force->numeric(arg[3]); double cut_one = cut_global; if (narg == 5) cut_one = force->numeric(arg[4]); if (cut_one <= 0.0) error->all(FLERR,"Incorrect args for pair coefficients"); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut[i][j] = cut_one; setflag[i][j] = 1; count++; } } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairLJShiftedForce::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] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); double ratio = sigma[i][j] / cut[i][j]; foffset[i][j] = 4.0 * epsilon[i][j] * (12.0 * pow(ratio,12.0) - 6.0 * pow(ratio,6.0)); offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); cut[j][i] = cut[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]; foffset[j][i] = foffset[i][j]; offset[j][i] = offset[i][j]; return cut[i][j]; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJShiftedForce::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { fwrite(&epsilon[i][j],sizeof(double),1,fp); fwrite(&sigma[i][j],sizeof(double),1,fp); fwrite(&cut[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJShiftedForce::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(&cut[i][j],sizeof(double),1,fp); } MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJShiftedForce::write_restart_settings(FILE *fp) { fwrite(&cut_global,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJShiftedForce::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { fread(&cut_global,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); } MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- */ double PairLJShiftedForce::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { double r2inv,r6inv,forcelj,philj,r,t; r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; r = sqrt(rsq); t = r/cut[itype][jtype]; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) - t*foffset[itype][jtype]; fforce = factor_lj*forcelj*r2inv; philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) + (t-1.0)*foffset[itype][jtype] - offset[itype][jtype]; return factor_lj*philj; } diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 2f0c735ef..49684020c 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -1,393 +1,389 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "math.h" #include "string.h" #include "stdlib.h" #include "fix_adapt.h" #include "atom.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}; /* ---------------------------------------------------------------------- */ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (narg < 5) error->all(FLERR,"Illegal fix adapt command"); nevery = atoi(arg[3]); if (nevery < 0) error->all(FLERR,"Illegal fix adapt 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 command"); nadapt++; iarg += 6; } else if (strcmp(arg[iarg],"kspace") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); nadapt++; iarg += 2; } else if (strcmp(arg[iarg],"atom") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command"); nadapt++; iarg += 3; } else break; } if (nadapt == 0) error->all(FLERR,"Illegal fix adapt 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 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 command"); nadapt++; iarg += 6; } else if (strcmp(arg[iarg],"kspace") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt 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 command"); nadapt++; iarg += 2; } else if (strcmp(arg[iarg],"atom") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt 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 command"); if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { int n = strlen(&arg[iarg+2][2]) + 1; adapt[nadapt].var = new char[n]; strcpy(adapt[nadapt].var,&arg[iarg+2][2]); } else error->all(FLERR,"Illegal fix adapt command"); nadapt++; iarg += 3; } else break; } // optional keywords resetflag = 0; scaleflag = 0; while (iarg < narg) { if (strcmp(arg[iarg],"reset") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt 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 command"); iarg += 2; } else if (strcmp(arg[iarg],"scale") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt 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 command"); iarg += 2; } else error->all(FLERR,"Illegal fix adapt 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"); } /* ---------------------------------------------------------------------- */ FixAdapt::~FixAdapt() { 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 FixAdapt::setmask() { int mask = 0; mask |= PRE_FORCE; mask |= POST_RUN; return mask; } /* ---------------------------------------------------------------------- */ void FixAdapt::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 does not exist"); if (!input->variable->equalstyle(ad->ivar)) error->all(FLERR,"Variable for fix adapt 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 pair style does not exist"); void *ptr = pair->extract(ad->pparam,ad->pdim); if (ptr == NULL) error->all(FLERR,"Fix adapt 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 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 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 requires atom attribute diameter"); } if (ad->aparam == CHARGE) { if (!atom->q_flag) error->all(FLERR,"Fix adapt 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]; } } } /* ---------------------------------------------------------------------- */ void FixAdapt::setup_pre_force(int vflag) { change_settings(); } /* ---------------------------------------------------------------------- */ void FixAdapt::pre_force(int vflag) { if (nevery == 0) return; if (update->ntimestep % nevery) return; change_settings(); } /* ---------------------------------------------------------------------- */ void FixAdapt::post_run() { if (resetflag) restore_settings(); } /* ---------------------------------------------------------------------- change pair,kspace,atom parameters based on variable evaluation ------------------------------------------------------------------------- */ void FixAdapt::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; } // 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; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; if (mflag == 0) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) radius[i] = 0.5*value; } else { for (i = 0; i < nlocal; i++) 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) { double *q = atom->q; int *mask = atom->mask; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { -// cout<<"Fix adapt called now!\n"; -// cout<<" Charge of atom: "<<i<<" is: "<<value; - q[i] = value; - } + if (mask[i] & groupbit) q[i] = value; } } } 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 FixAdapt::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(); }