diff --git a/src/CORESHELL/compute_temp_cs.cpp b/src/CORESHELL/compute_temp_cs.cpp index 6d0aa2514a..7512d48afd 100644 --- a/src/CORESHELL/compute_temp_cs.cpp +++ b/src/CORESHELL/compute_temp_cs.cpp @@ -54,9 +54,6 @@ ComputeTempCS::ComputeTempCS(LAMMPS *lmp, int narg, char **arg) : tempbias = 1; extarray = 0; - int *mask = atom->mask; - int nlocal = atom->nlocal; - // find and define groupbits for core and shell groups cgroup = group->find(arg[3]); @@ -226,7 +223,6 @@ void ComputeTempCS::dof_compute() double ComputeTempCS::compute_scalar() { - int i; double vthermal[3]; invoked_scalar = update->ntimestep; @@ -280,7 +276,6 @@ void ComputeTempCS::compute_vector() double **v = atom->v; int *mask = atom->mask; - tagint *molecule = atom->molecule; int *type = atom->type; double *mass = atom->mass; double *rmass = atom->rmass; @@ -340,7 +335,7 @@ void ComputeTempCS::vcm_pairs() double *partner = fix->vstore; tagint partnerID; - for (int i = 0; i < nlocal; i++) { + for (i = 0; i < nlocal; i++) { if ((mask[i] & groupbit) && (mask[i] & groupbit_c || mask[i] & groupbit_s)) { if (rmass) massone = rmass[i]; diff --git a/src/USER-FEP/compute_fep.h b/src/USER-FEP/compute_fep.h index 7c60a68f71..51fcbd5ea6 100644 --- a/src/USER-FEP/compute_fep.h +++ b/src/USER-FEP/compute_fep.h @@ -42,7 +42,6 @@ class ComputeFEP : public Compute { int tailflag, volumeflag; int fepinitflag; int eflag, vflag; - int sys_qsum_update_flag; double temp_fep; int nmax; diff --git a/src/USER-MOLFILE/dump_molfile.cpp b/src/USER-MOLFILE/dump_molfile.cpp index b0e4d54f61..add7bf695b 100644 --- a/src/USER-MOLFILE/dump_molfile.cpp +++ b/src/USER-MOLFILE/dump_molfile.cpp @@ -81,7 +81,7 @@ DumpMolfile::DumpMolfile(LAMMPS *lmp, int narg, char **arg) // allocate global array for atom coords bigint n = group->count(igroup); - if (n > MAXSMALLINT/sizeof(float)) + if (n > static_cast(MAXSMALLINT/3/sizeof(float))) error->all(FLERR,"Too many atoms for dump molfile"); if (n < 1) error->all(FLERR,"Not enough atoms for dump molfile"); diff --git a/src/USER-OMP/Install.sh b/src/USER-OMP/Install.sh index 547fa074e6..85b44f1bee 100644 --- a/src/USER-OMP/Install.sh +++ b/src/USER-OMP/Install.sh @@ -51,11 +51,14 @@ if (test $mode = 1) then sed -i -e 's|^PKG_INC =[ \t]*|&-DLMP_USER_OMP |' ../Makefile.package fi - # need to delete a bunch of depenency files because they indirectly depend on user_cuda.h + # need to delete a bunch of dependency files because they + # indirectly depend on user_cuda.h + for f in finish.d modify_cuda.d do \ rm -f ../Obj_*/$f done + # force rebuild of files with LMP_USER_OMP switch touch ../accelerator_omp.h @@ -66,6 +69,14 @@ elif (test $mode = 0) then sed -i -e 's/[^ \t]*OMP[^ \t]* //' ../Makefile.package fi + # need to delete a bunch of dependency files because they + # indirectly depend on user_cuda.h + + for f in finish.d modify_cuda.d + do \ + rm -f ../Obj_*/$f + done + # force rebuild of files with LMP_USER_OMP switch touch ../accelerator_omp.h diff --git a/src/XTC/dump_xtc.cpp b/src/XTC/dump_xtc.cpp index 5f7bee27db..7b8afbd3a9 100644 --- a/src/XTC/dump_xtc.cpp +++ b/src/XTC/dump_xtc.cpp @@ -70,7 +70,7 @@ DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) // allocate global array for atom coords bigint n = group->count(igroup); - if (n > MAXSMALLINT/3/sizeof(float)) + if (n > static_cast(MAXSMALLINT/3/sizeof(float))) error->all(FLERR,"Too many atoms for dump xtc"); natoms = static_cast (n); diff --git a/src/domain.cpp b/src/domain.cpp index 005678c952..da93f50720 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -451,9 +451,13 @@ void Domain::reset_box() set_global_box(); set_local_box(); - // if shrink-wrapped & kspace is defined (i.e. using MSM) call setup() + // if shrink-wrapped & kspace is defined (i.e. using MSM), call setup() + // also call init() (to test for compatibility) ? - if (nonperiodic == 2 && force->kspace) force->kspace->setup(); + if (nonperiodic == 2 && force->kspace) { + //force->kspace->init(); + force->kspace->setup(); + } // if shrink-wrapped & triclinic, re-convert to lamda coords for new box // re-invoke pbc() b/c x2lamda result can be outside [0,1] due to roundoff diff --git a/src/dump_dcd.cpp b/src/dump_dcd.cpp index a4a6e4f704..c679df62b9 100644 --- a/src/dump_dcd.cpp +++ b/src/dump_dcd.cpp @@ -68,7 +68,7 @@ DumpDCD::DumpDCD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) // allocate global array for atom coords bigint n = group->count(igroup); - if (n > MAXSMALLINT/sizeof(float)) + if (n > static_cast(MAXSMALLINT/3/sizeof(float))) error->all(FLERR,"Too many atoms for dump dcd"); natoms = static_cast (n); diff --git a/src/fix_ave_chunk.cpp b/src/fix_ave_chunk.cpp index 2eb3f9dce1..7a7935c5bb 100644 --- a/src/fix_ave_chunk.cpp +++ b/src/fix_ave_chunk.cpp @@ -574,7 +574,6 @@ void FixAveChunk::end_of_step() // sum within each chunk, only include atoms in fix group // compute/fix/variable may invoke computes so wrap with clear/add - double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; diff --git a/src/fix_temp_csld.cpp b/src/fix_temp_csld.cpp new file mode 100644 index 0000000000..e394f06e58 --- /dev/null +++ b/src/fix_temp_csld.cpp @@ -0,0 +1,307 @@ +/* ---------------------------------------------------------------------- + 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: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "stdlib.h" +#include "math.h" +#include "fix_temp_csld.h" +#include "atom.h" +#include "force.h" +#include "memory.h" +#include "comm.h" +#include "input.h" +#include "variable.h" +#include "group.h" +#include "update.h" +#include "modify.h" +#include "compute.h" +#include "random_mars.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{NOBIAS,BIAS}; +enum{CONSTANT,EQUAL}; + +/* ---------------------------------------------------------------------- */ + +FixTempCSLD::FixTempCSLD(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg != 7) error->all(FLERR,"Illegal fix temp/csld command"); + + // CSLD thermostat should be applied every step + + nevery = 1; + scalar_flag = 1; + global_freq = nevery; + dynamic_group_allow = 1; + extscalar = 1; + + tstr = NULL; + if (strstr(arg[3],"v_") == arg[3]) { + int n = strlen(&arg[3][2]) + 1; + tstr = new char[n]; + strcpy(tstr,&arg[3][2]); + tstyle = EQUAL; + } else { + t_start = force->numeric(FLERR,arg[3]); + t_target = t_start; + tstyle = CONSTANT; + } + + t_stop = force->numeric(FLERR,arg[4]); + t_period = force->numeric(FLERR,arg[5]); + int seed = force->inumeric(FLERR,arg[6]); + + // error checks + + if (t_period <= 0.0) error->all(FLERR,"Illegal fix temp/csld command"); + if (seed <= 0) error->all(FLERR,"Illegal fix temp/csld command"); + + random = new RanMars(lmp,seed + comm->me); + + // create a new compute temp style + // id = fix-ID + temp, compute group = fix group + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = group->names[igroup]; + newarg[2] = (char *) "temp"; + modify->add_compute(3,newarg); + delete [] newarg; + tflag = 1; + + vhold = NULL; + nmax = -1; + energy = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +FixTempCSLD::~FixTempCSLD() +{ + delete [] tstr; + + // delete temperature if fix created it + + if (tflag) modify->delete_compute(id_temp); + delete [] id_temp; + + delete random; + memory->destroy(vhold); + vhold = NULL; + nmax = -1; +} + +/* ---------------------------------------------------------------------- */ + +int FixTempCSLD::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + mask |= THERMO_ENERGY; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixTempCSLD::init() +{ + + // we cannot handle constraints via rattle or shake correctly. + + int has_shake = 0; + for (int i = 0; i < modify->nfix; i++) + if ((strcmp(modify->fix[i]->style,"shake") == 0) + || (strcmp(modify->fix[i]->style,"rattle") == 0)) ++has_shake; + + if (has_shake > 0) + error->all(FLERR,"Fix temp/csld is not compatible with fix rattle or fix shake"); + + // check variable + + if (tstr) { + tvar = input->variable->find(tstr); + if (tvar < 0) + error->all(FLERR,"Variable name for fix temp/csld does not exist"); + if (input->variable->equalstyle(tvar)) tstyle = EQUAL; + else error->all(FLERR,"Variable for fix temp/csld is invalid style"); + } + + int icompute = modify->find_compute(id_temp); + if (icompute < 0) + error->all(FLERR,"Temperature ID for fix temp/csld does not exist"); + temperature = modify->compute[icompute]; + + if (temperature->tempbias) which = BIAS; + else which = NOBIAS; +} + +/* ---------------------------------------------------------------------- */ + +void FixTempCSLD::end_of_step() +{ + + // set current t_target + // if variable temp, evaluate variable, wrap with clear/add + + double delta = update->ntimestep - update->beginstep; + + if (delta != 0.0) delta /= update->endstep - update->beginstep; + if (tstyle == CONSTANT) + t_target = t_start + delta * (t_stop-t_start); + else { + modify->clearstep_compute(); + t_target = input->variable->compute_equal(tvar); + if (t_target < 0.0) + error->one(FLERR, + "Fix temp/csld variable returned negative temperature"); + modify->addstep_compute(update->ntimestep + nevery); + } + + double t_current = temperature->compute_scalar(); + double ekin_old = t_current * 0.5 * temperature->dof * force->boltz; + + double * const * const v = atom->v; + const int * const mask = atom->mask; + const int * const type = atom->type; + const int nlocal = atom->nlocal; + + // adjust holding space, if needed and copy existing velocities + + if (nmax < atom->nlocal) { + nmax = atom->nlocal + 1; + memory->destroy(vhold); + memory->create(vhold,nmax,3,"csld:vhold"); + } + + // The CSLD thermostat is a linear combination of old and new velocities, + // where the new ones are randomly chosen from a gaussian distribution. + // see Bussi and Parrinello, Phys. Rev. E (2007). + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + double m; + if (atom->rmass_flag) m = atom->rmass[i]; + else m = atom->mass[type[i]]; + + const double factor = 1.0/sqrt(m); + const double vx = random->gaussian() * factor; + vhold[i][0] = v[i][0]; + v[i][0] = vx; + const double vy = random->gaussian() * factor; + vhold[i][1] = v[i][1]; + v[i][1] = vy; + const double vz = random->gaussian() * factor; + vhold[i][2] = v[i][2]; + v[i][2] = vz; + } + } + + // mixing factors + const double c1 = exp(-update->dt/t_period); + const double c2 = sqrt((1.0-c1*c1)*t_target/temperature->compute_scalar()); + + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + v[i][0] = vhold[i][0]*c1 + v[i][0]*c2; + v[i][1] = vhold[i][1]*c1 + v[i][1]*c2; + v[i][2] = vhold[i][2]*c1 + v[i][2]*c2; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,vhold[i]); + v[i][0] = vhold[i][0]*c1 + v[i][0]*c2; + v[i][1] = vhold[i][1]*c1 + v[i][1]*c2; + v[i][2] = vhold[i][2]*c1 + v[i][2]*c2; + temperature->restore_bias(i,v[i]); + } + } + } + + // tally the kinetic energy transferred between heat bath and system + + t_current = temperature->compute_scalar(); + energy += ekin_old - t_current * 0.5 * temperature->dof * force->boltz; +} + +/* ---------------------------------------------------------------------- */ + +int FixTempCSLD::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"temp") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (tflag) { + modify->delete_compute(id_temp); + tflag = 0; + } + delete [] id_temp; + int n = strlen(arg[1]) + 1; + id_temp = new char[n]; + strcpy(id_temp,arg[1]); + + int icompute = modify->find_compute(id_temp); + if (icompute < 0) + error->all(FLERR,"Could not find fix_modify temperature ID"); + temperature = modify->compute[icompute]; + + if (temperature->tempflag == 0) + error->all(FLERR, + "Fix_modify temperature ID does not compute temperature"); + if (temperature->igroup != igroup && comm->me == 0) + error->warning(FLERR,"Group for fix_modify temp != fix group"); + return 2; + } + return 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixTempCSLD::reset_target(double t_new) +{ + t_target = t_start = t_stop = t_new; +} + +/* ---------------------------------------------------------------------- */ + +double FixTempCSLD::compute_scalar() +{ + return energy; +} + +/* ---------------------------------------------------------------------- + extract thermostat properties +------------------------------------------------------------------------- */ + +void *FixTempCSLD::extract(const char *str, int &dim) +{ + dim=0; + if (strcmp(str,"t_target") == 0) { + return &t_target; + } + return NULL; +} diff --git a/src/fix_temp_csld.h b/src/fix_temp_csld.h new file mode 100644 index 0000000000..3c4ddeefb4 --- /dev/null +++ b/src/fix_temp_csld.h @@ -0,0 +1,101 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(temp/csld,FixTempCSLD) + +#else + +#ifndef LMP_FIX_TEMP_CSLD_H +#define LMP_FIX_TEMP_CSLD_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixTempCSLD : public Fix { + public: + FixTempCSLD(class LAMMPS *, int, char **); + ~FixTempCSLD(); + int setmask(); + void init(); + void end_of_step(); + int modify_param(int, char **); + void reset_target(double); + virtual double compute_scalar(); + virtual void *extract(const char *, int &); + + private: + double t_start,t_stop,t_period,t_target; + double **vhold; + double energy; + int nmax,which; + int tstyle,tvar; + char *tstr; + + char *id_temp; + class Compute *temperature; + int tflag; + + class RanMars *random; +}; + +} + +#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: Fix temp/csld is not compatible with fix rattle or fix shake + +These two commands cannot currently be used together with fix temp/csld. + +E: Variable name for fix temp/csld does not exist + +Self-explanatory. + +E: Variable for fix temp/csld is invalid style + +Only equal-style variables can be used. + +E: Temperature ID for fix temp/csld does not exist + +Self-explanatory. + +E: Fix temp/csld variable returned negative temperature + +Self-explanatory. + +E: Could not find fix_modify temperature ID + +The compute ID for computing temperature does not exist. + +E: Fix_modify temperature ID does not compute temperature + +The compute ID assigned to the fix must compute temperature. + +W: Group for fix_modify temp != fix group + +The fix_modify command is specifying a temperature computation that +computes a temperature on a different group of atoms than the fix +itself operates on. This is probably not what you want to do. + +*/ diff --git a/src/fix_temp_csvr.cpp b/src/fix_temp_csvr.cpp index a8c0d9b508..09ffb951dd 100644 --- a/src/fix_temp_csvr.cpp +++ b/src/fix_temp_csvr.cpp @@ -12,7 +12,8 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Axel Kohlmeyer (ICTP, Italy) + Contributing author: Axel Kohlmeyer (Temple U) + Based on code by Paulo Raiteri (Curtin U) and Giovanni Bussi (SISSA) ------------------------------------------------------------------------- */ #include "string.h" @@ -38,6 +39,78 @@ using namespace FixConst; enum{NOBIAS,BIAS}; enum{CONSTANT,EQUAL}; +double FixTempCSVR::gamdev(const int ia) +{ + int j; + double am,e,s,v1,v2,x,y; + + if (ia < 1) return 0.0; + if (ia < 6) { + x=1.0; + for (j=1; j<=ia; j++) + x *= random->uniform(); + x = -log(x); + } else { + restart: + do { + do { + do { + v1 = random->uniform(); + v2 = 2.0*random->uniform() - 1.0; + } while (v1*v1 + v2*v2 > 1.0); + + y=v2/v1; + am=ia-1; + s=sqrt(2.0*am+1.0); + x=s*y+am; + } while (x <= 0.0); + + if (am*log(x/am)-s*y < -700 || v1<0.00001) { + goto restart; + } + + e=(1.0+y*y)*exp(am*log(x/am)-s*y); + } while (random->uniform() > e); + } + return x; +} + +/* ------------------------------------------------------------------- + returns the sum of n independent gaussian noises squared + (i.e. equivalent to summing the square of the return values of nn + calls to gasdev) +---------------------------------------------------------------------- */ +double FixTempCSVR::sumnoises(int nn) { + if (nn == 0) { + return 0.0; + } else if (nn == 1) { + const double rr = random->gaussian(); + return rr*rr; + } else if (nn % 2 == 0) { + return 2.0 * gamdev(nn / 2); + } else { + const double rr = random->gaussian(); + return 2.0 * gamdev((nn-1) / 2) + rr*rr; + } + return 0.0; +} + +/* ------------------------------------------------------------------- + returns the scaling factor for velocities to thermalize + the system so it samples the canonical ensemble +---------------------------------------------------------------------- */ + +double FixTempCSVR::resamplekin(double ekin_old, double ekin_new){ + const double tdof = temperature->dof; + const double c1 = exp(-update->dt/t_period); + const double c2 = (1.0-c1)*ekin_new/ekin_old/tdof; + const double r1 = random->gaussian(); + const double r2 = sumnoises(tdof - 1); + + const double scale = c1 + c2*(r1*r1+r2) + 2.0*r1*sqrt(c1*c2); + return sqrt(scale); +} + /* ---------------------------------------------------------------------- */ FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) : @@ -48,8 +121,10 @@ FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) : // CSVR thermostat should be applied every step nevery = 1; + scalar_flag = 1; global_freq = nevery; dynamic_group_allow = 1; + extscalar = 1; tstr = NULL; if (strstr(arg[3],"v_") == arg[3]) { @@ -72,6 +147,8 @@ FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) : if (t_period <= 0.0) error->all(FLERR,"Illegal fix temp/csvr command"); if (seed <= 0) error->all(FLERR,"Illegal fix temp/csvr command"); + random = new RanMars(lmp,seed + comm->me); + // create a new compute temp style // id = fix-ID + temp, compute group = fix group @@ -88,9 +165,6 @@ FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) : delete [] newarg; tflag = 1; - random = new RanMars(lmp,seed + comm->me); - - vhold = NULL; nmax = -1; energy = 0.0; } @@ -107,8 +181,6 @@ FixTempCSVR::~FixTempCSVR() delete [] id_temp; delete random; - memory->destroy(vhold); - vhold = NULL; nmax = -1; } @@ -118,6 +190,7 @@ int FixTempCSVR::setmask() { int mask = 0; mask |= END_OF_STEP; + mask |= THERMO_ENERGY; return mask; } @@ -125,15 +198,7 @@ int FixTempCSVR::setmask() void FixTempCSVR::init() { - // we cannot handle shake correctly at the moment - int has_shake = 0; - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"shake") == 0) ++has_shake; - - if (has_shake > 0) - error->all(FLERR,"Fix temp/csvr is not compatible with fix shake"); - // check variable if (tstr) { @@ -175,65 +240,45 @@ void FixTempCSVR::end_of_step() modify->addstep_compute(update->ntimestep + nevery); } + const double t_current = temperature->compute_scalar(); + const double efactor = 0.5 * temperature->dof * force->boltz; + const double ekin_old = t_current * efactor; + const double ekin_new = t_target * efactor; + + // compute velocity scaling factor on root node and broadcast + double lamda; + if (comm->me == 0) { + lamda = resamplekin(ekin_old, ekin_new); + } + MPI_Bcast(&lamda,1,MPI_DOUBLE,0,world); + double * const * const v = atom->v; const int * const mask = atom->mask; - const int * const type = atom->type; const int nlocal = atom->nlocal; - // adjust holding space, if needed and copy existing velocities - - if (nmax < atom->nlocal) { - nmax = atom->nlocal + 1; - memory->destroy(vhold); - memory->create(vhold,nmax,3,"csvr:vhold"); - } - - // The CSVR thermostat is a linear combination of old and new velocities, - // where the new ones are randomly chosen from a gaussian distribution. - // see Bussi and Parrinello, Phys. Rev. E (2007). - - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - double m; - if (atom->rmass_flag) m = atom->rmass[i]; - else m = atom->mass[type[i]]; - - const double factor = 1.0/sqrt(m); - const double vx = random->gaussian() * factor; - vhold[i][0] = v[i][0]; - v[i][0] = vx; - const double vy = random->gaussian() * factor; - vhold[i][1] = v[i][1]; - v[i][1] = vy; - const double vz = random->gaussian() * factor; - vhold[i][2] = v[i][2]; - v[i][2] = vz; - } - } - - // mixing factors - const double c1 = exp(-update->dt/t_period); - const double c2 = sqrt((1.0-c1*c1)*t_target/temperature->compute_scalar()); - if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - v[i][0] = vhold[i][0]*c1 + v[i][0]*c2; - v[i][1] = vhold[i][1]*c1 + v[i][1]*c2; - v[i][2] = vhold[i][2]*c1 + v[i][2]*c2; + v[i][0] *= lamda; + v[i][1] *= lamda; + v[i][2] *= lamda; } } } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - temperature->remove_bias(i,vhold[i]); - v[i][0] = vhold[i][0]*c1 + v[i][0]*c2; - v[i][1] = vhold[i][1]*c1 + v[i][1]*c2; - v[i][2] = vhold[i][2]*c1 + v[i][2]*c2; + temperature->remove_bias(i,v[i]); + v[i][0] *= lamda; + v[i][1] *= lamda; + v[i][2] *= lamda; temperature->restore_bias(i,v[i]); } } } + + // tally the kinetic energy transferred between heat bath and system + + energy += ekin_old * (1.0 - lamda*lamda); } /* ---------------------------------------------------------------------- */ @@ -273,6 +318,13 @@ void FixTempCSVR::reset_target(double t_new) t_target = t_start = t_stop = t_new; } +/* ---------------------------------------------------------------------- */ + +double FixTempCSVR::compute_scalar() +{ + return energy; +} + /* ---------------------------------------------------------------------- extract thermostat properties ------------------------------------------------------------------------- */ diff --git a/src/fix_temp_csvr.h b/src/fix_temp_csvr.h index 8fcdddb404..48e7cbce8e 100644 --- a/src/fix_temp_csvr.h +++ b/src/fix_temp_csvr.h @@ -33,11 +33,11 @@ class FixTempCSVR : public Fix { void end_of_step(); int modify_param(int, char **); void reset_target(double); + virtual double compute_scalar(); virtual void *extract(const char *, int &); private: double t_start,t_stop,t_period,t_target; - double **vhold; double energy; int nmax,which; int tstyle,tvar; @@ -48,6 +48,11 @@ class FixTempCSVR : public Fix { int tflag; class RanMars *random; + + private: + double resamplekin(double, double); + double sumnoises(int); + double gamdev(int); }; } @@ -63,9 +68,6 @@ 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: Fix temp/csvr is not compatible with fix shake - -These two commands cannot currently be used toghether. E: Variable name for fix temp/csvr does not exist diff --git a/src/force.cpp b/src/force.cpp index 71eb0add9d..ac6f455321 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -677,21 +677,21 @@ void Force::set_special(int narg, char **arg) iarg += 1; } else if (strcmp(arg[iarg],"lj/coul") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal special_bonds command"); - special_lj[1] = special_coul[1] = atof(arg[iarg+1]); - special_lj[2] = special_coul[2] = atof(arg[iarg+2]); - special_lj[3] = special_coul[3] = atof(arg[iarg+3]); + special_lj[1] = special_coul[1] = numeric(FLERR,arg[iarg+1]); + special_lj[2] = special_coul[2] = numeric(FLERR,arg[iarg+2]); + special_lj[3] = special_coul[3] = numeric(FLERR,arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg],"lj") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal special_bonds command"); - special_lj[1] = atof(arg[iarg+1]); - special_lj[2] = atof(arg[iarg+2]); - special_lj[3] = atof(arg[iarg+3]); + special_lj[1] = numeric(FLERR,arg[iarg+1]); + special_lj[2] = numeric(FLERR,arg[iarg+2]); + special_lj[3] = numeric(FLERR,arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg],"coul") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal special_bonds command"); - special_coul[1] = atof(arg[iarg+1]); - special_coul[2] = atof(arg[iarg+2]); - special_coul[3] = atof(arg[iarg+3]); + special_coul[1] = numeric(FLERR,arg[iarg+1]); + special_coul[2] = numeric(FLERR,arg[iarg+2]); + special_coul[3] = numeric(FLERR,arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg],"angle") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal special_bonds command"); diff --git a/src/variable.cpp b/src/variable.cpp index f9828a4a12..df14be27a8 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -2308,8 +2308,8 @@ double Variable::collapse_tree(Tree *tree) else if (update->ntimestep < ivalue2) { int offset = update->ntimestep - ivalue1; tree->value = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3; - if (tree->value > ivalue2) tree->value = 9.0e18; - } else tree->value = 9.0e18; + if (tree->value > ivalue2) tree->value = MAXBIGINT; + } else tree->value = MAXBIGINT; return tree->value; } @@ -2345,10 +2345,10 @@ double Variable::collapse_tree(Tree *tree) if (istep > ivalue5) { int offset = ivalue5 - ivalue1; istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3; - if (istep > ivalue2) istep = 9.0e18; + if (istep > ivalue2) istep = MAXBIGINT; } } - } else istep = 9.0e18; + } else istep = MAXBIGINT; tree->value = istep; return tree->value; } @@ -2633,8 +2633,8 @@ double Variable::eval_tree(Tree *tree, int i) else if (update->ntimestep < ivalue2) { int offset = update->ntimestep - ivalue1; arg = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3; - if (arg > ivalue2) arg = 9.0e18; - } else arg = 9.0e18; + if (arg > ivalue2) arg = MAXBIGINT; + } else arg = MAXBIGINT; return arg; } @@ -2665,10 +2665,10 @@ double Variable::eval_tree(Tree *tree, int i) if (istep > ivalue5) { int offset = ivalue5 - ivalue1; istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3; - if (istep > ivalue2) istep = 9.0e18; + if (istep > ivalue2) istep = MAXBIGINT; } } - } else istep = 9.0e18; + } else istep = MAXBIGINT; arg = istep; return arg; } @@ -3141,8 +3141,8 @@ int Variable::math_function(char *word, char *contents, Tree **tree, else if (update->ntimestep < ivalue2) { int offset = update->ntimestep - ivalue1; value = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3; - if (value > ivalue2) value = 9.0e18; - } else value = 9.0e18; + if (value > ivalue2) value = MAXBIGINT; + } else value = MAXBIGINT; argstack[nargstack++] = value; } @@ -3176,10 +3176,10 @@ int Variable::math_function(char *word, char *contents, Tree **tree, if (istep > ivalue5) { int offset = ivalue5 - ivalue1; istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3; - if (istep > ivalue2) istep = 9.0e18; + if (istep > ivalue2) istep = MAXBIGINT; } } - } else istep = 9.0e18; + } else istep = MAXBIGINT; double value = istep; argstack[nargstack++] = value; }