jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13320 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2015-03-27 20:39:27 +00:00
parent bc06bcf854
commit 9435c9dcb1
14 changed files with 566 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
src/CORESHELL/compute_temp_cs.cpp
View File

@ -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];

1
src/USER-FEP/compute_fep.h
View File

@ -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;

2
src/USER-MOLFILE/dump_molfile.cpp
View File

@ -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<bigint>(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");

13
src/USER-OMP/Install.sh
View File

@ -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

2
src/XTC/dump_xtc.cpp
View File

@ -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<bigint>(MAXSMALLINT/3/sizeof(float)))
error->all(FLERR,"Too many atoms for dump xtc");
natoms = static_cast<int> (n);

8
src/domain.cpp
View File

@ -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

2
src/dump_dcd.cpp
View File

@ -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<bigint>(MAXSMALLINT/3/sizeof(float)))
error->all(FLERR,"Too many atoms for dump dcd");
natoms = static_cast<int> (n);

1
src/fix_ave_chunk.cpp
View File

@ -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;

307
src/fix_temp_csld.cpp Normal file
View File

@ -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;
}

101
src/fix_temp_csld.h Normal file
View File

@ -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.
*/

166
src/fix_temp_csvr.cpp
View File

@ -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
------------------------------------------------------------------------- */

10
src/fix_temp_csvr.h
View File

@ -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

18
src/force.cpp
View File

@ -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");

24
src/variable.cpp
View File

@ -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;
}