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

This commit is contained in:
sjplimp
2013-11-06 19:00:54 +00:00
parent aa9502ec5a
commit faff0eb70b
4 changed files with 722 additions and 704 deletions

View File

@ -1,197 +1,197 @@
/* ------------------------------------------------------------------------- /* -------------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------- /* -------------------------------------------------------------------------
Contributing authors: Contributing authors:
Rodrigo Freitas (Unicamp/Brazil) - rodrigohb@gmail.com Rodrigo Freitas (Unicamp/Brazil) - rodrigohb@gmail.com
Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "stdlib.h" #include "stdlib.h"
#include "string.h" #include "string.h"
#include "math.h" #include "math.h"
#include "fix_ti_rs.h" #include "fix_ti_rs.h"
#include "atom.h" #include "atom.h"
#include "update.h" #include "update.h"
#include "respa.h" #include "respa.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Class constructor initialize all variables. // Class constructor initialize all variables.
FixTIRS::FixTIRS(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) FixTIRS::FixTIRS(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{ {
// Checking the input information. // Checking the input information.
if (narg < 7 || narg > 9) error->all(FLERR,"Illegal fix ti/rs command"); if (narg < 7 || narg > 9) error->all(FLERR,"Illegal fix ti/rs command");
// Fix flags. // Fix flags.
vector_flag = 1; vector_flag = 1;
size_vector = 2; size_vector = 2;
global_freq = 1; global_freq = 1;
extvector = 1; extvector = 1;
// Time variables. // Time variables.
t_switch = atof(arg[5]); t_switch = atoi(arg[5]);
t_equil = atof(arg[6]); t_equil = atoi(arg[6]);
t0 = update->ntimestep; t0 = update->ntimestep;
if (t_switch < 0.0) error->all(FLERR,"Illegal fix ti/rs command"); if (t_switch < 0.0) error->all(FLERR,"Illegal fix ti/rs command");
if (t_equil < 0.0) error->all(FLERR,"Illegal fix ti/rs command"); if (t_equil < 0.0) error->all(FLERR,"Illegal fix ti/rs command");
// Coupling parameter limits and initialization. // Coupling parameter limits and initialization.
l_initial = atof(arg[3]); l_initial = atof(arg[3]);
l_final = atof(arg[4]); l_final = atof(arg[4]);
sf = 1; sf = 1;
if (narg > 7) { if (narg > 7) {
if (strcmp(arg[7], "function") == 0) sf = atoi(arg[8]); if (strcmp(arg[7], "function") == 0) sf = atoi(arg[8]);
else error->all(FLERR,"Illegal fix ti/rs switching function"); else error->all(FLERR,"Illegal fix ti/rs switching function");
if ((sf<1) || (sf>3)) if ((sf<1) || (sf>3))
error->all(FLERR,"Illegal fix ti/rs switching function"); error->all(FLERR,"Illegal fix ti/rs switching function");
} }
lambda = switch_func(0); lambda = switch_func(0);
dlambda = dswitch_func(0); dlambda = dswitch_func(0);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixTIRS::~FixTIRS() FixTIRS::~FixTIRS()
{ {
// unregister callbacks to this fix from Atom class // unregister callbacks to this fix from Atom class
atom->delete_callback(id,0); atom->delete_callback(id,0);
atom->delete_callback(id,1); atom->delete_callback(id,1);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
int FixTIRS::setmask() int FixTIRS::setmask()
{ {
int mask = 0; int mask = 0;
mask |= INITIAL_INTEGRATE; mask |= INITIAL_INTEGRATE;
mask |= POST_FORCE; mask |= POST_FORCE;
return mask; return mask;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTIRS::init() void FixTIRS::init()
{ {
if (strstr(update->integrate_style,"respa")) if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels; nlevels_respa = ((Respa *) update->integrate)->nlevels;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTIRS::setup(int vflag) void FixTIRS::setup(int vflag)
{ {
if (strstr(update->integrate_style,"verlet")) if (strstr(update->integrate_style,"verlet"))
post_force(vflag); post_force(vflag);
else { else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0); post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
} }
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTIRS::min_setup(int vflag) void FixTIRS::min_setup(int vflag)
{ {
post_force(vflag); post_force(vflag);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTIRS::post_force(int vflag) void FixTIRS::post_force(int vflag)
{ {
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double **f = atom->f; double **f = atom->f;
// Scaling forces. // Scaling forces.
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
f[i][0] = lambda * f[i][0]; f[i][0] = lambda * f[i][0];
f[i][1] = lambda * f[i][1]; f[i][1] = lambda * f[i][1];
f[i][2] = lambda * f[i][2]; f[i][2] = lambda * f[i][2];
} }
} }
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTIRS::post_force_respa(int vflag, int ilevel, int iloop) void FixTIRS::post_force_respa(int vflag, int ilevel, int iloop)
{ {
if (ilevel == nlevels_respa-1) post_force(vflag); if (ilevel == nlevels_respa-1) post_force(vflag);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTIRS::min_post_force(int vflag) void FixTIRS::min_post_force(int vflag)
{ {
post_force(vflag); post_force(vflag);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTIRS::initial_integrate(int vflag) void FixTIRS::initial_integrate(int vflag)
{ {
// Update the coupling parameter value. // Update the coupling parameter value.
double t = update->ntimestep - (t0+t_equil); double t = update->ntimestep - (t0+t_equil);
if( (t >= 0) && (t <= t_switch) ) { if( (t >= 0) && (t <= t_switch) ) {
lambda = switch_func(t/t_switch); lambda = switch_func(t/t_switch);
dlambda = dswitch_func(t/t_switch); dlambda = dswitch_func(t/t_switch);
} }
if( (t >= t_equil+t_switch) && (t <= (t_equil+2*t_switch)) ) { if( (t >= t_equil+t_switch) && (t <= (t_equil+2*t_switch)) ) {
lambda = switch_func(1.0 - (t - t_switch - t_equil)/t_switch); lambda = switch_func(1.0 - (t - t_switch - t_equil)/t_switch);
dlambda = - dswitch_func(1.0 - (t - t_switch - t_equil)/t_switch); dlambda = - dswitch_func(1.0 - (t - t_switch - t_equil)/t_switch);
} }
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double FixTIRS::compute_vector(int n) double FixTIRS::compute_vector(int n)
{ {
linfo[0] = lambda; linfo[0] = lambda;
linfo[1] = dlambda; linfo[1] = dlambda;
return linfo[n]; return linfo[n];
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double FixTIRS::switch_func(double t) double FixTIRS::switch_func(double t)
{ {
if (sf == 2) return l_initial / (1 + t * (l_initial/l_final - 1)); if (sf == 2) return l_initial / (1 + t * (l_initial/l_final - 1));
if (sf == 3) return l_initial / (1 + log2(1+t) * (l_initial/l_final - 1)); if (sf == 3) return l_initial / (1 + log2(1+t) * (l_initial/l_final - 1));
// Default option is sf = 1. // Default option is sf = 1.
return l_initial + (l_final - l_initial) * t; return l_initial + (l_final - l_initial) * t;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double FixTIRS::dswitch_func(double t) double FixTIRS::dswitch_func(double t)
{ {
double aux = (1.0/l_initial - 1.0/l_final); double aux = (1.0/l_initial - 1.0/l_final);
if (sf == 2) return lambda * lambda * aux / t_switch; if (sf == 2) return lambda * lambda * aux / t_switch;
if (sf == 3) return lambda * lambda * aux / (t_switch * log(2) * (1 + t)); if (sf == 3) return lambda * lambda * aux / (t_switch * log(2) * (1 + t));
// Default option is sf = 1. // Default option is sf = 1.
return (l_final-l_initial)/t_switch; return (l_final-l_initial)/t_switch;
} }

View File

@ -1,60 +1,66 @@
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#ifdef FIX_CLASS /* ----------------------------------------------------------------------
Contributing authors:
FixStyle(ti/rs,FixTIRS) Rodrigo Freitas (Unicamp/Brazil) - rodrigohb@gmail.com
Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br
#else ------------------------------------------------------------------------- */
#ifndef LMP_FIX_TI_RS_H #ifdef FIX_CLASS
#define LMP_FIX_TI_RS_H
FixStyle(ti/rs,FixTIRS)
#include "fix.h"
#else
namespace LAMMPS_NS {
#ifndef LMP_FIX_TI_RS_H
class FixTIRS : public Fix { #define LMP_FIX_TI_RS_H
public:
FixTIRS(class LAMMPS *, int, char **); #include "fix.h"
~FixTIRS();
int setmask(); namespace LAMMPS_NS {
void init();
void setup(int); class FixTIRS : public Fix {
void min_setup(int); public:
void post_force(int); FixTIRS(class LAMMPS *, int, char **);
void post_force_respa(int, int, int); ~FixTIRS();
void min_post_force(int); int setmask();
virtual void initial_integrate(int); void init();
double compute_vector(int); void setup(int);
void min_setup(int);
private: void post_force(int);
double switch_func(double); void post_force_respa(int, int, int);
double dswitch_func(double); void min_post_force(int);
virtual void initial_integrate(int);
double lambda; // Coupling parameter. double compute_vector(int);
double dlambda; // Lambda variation with t.
double l_initial; // Lambda initial value. private:
double l_final; // Lambda final value. double switch_func(double);
double linfo[2]; // Current lambda status. double dswitch_func(double);
int t_switch; // Total switching steps.
int t_equil; // Equilibration time. double lambda; // Coupling parameter.
int t0; // Initial time. double dlambda; // Lambda variation with t.
int sf; // Switching function option. double l_initial; // Lambda initial value.
int nlevels_respa; double l_final; // Lambda final value.
}; double linfo[2]; // Current lambda status.
int t_switch; // Total switching steps.
} int t_equil; // Equilibration time.
int t0; // Initial time.
#endif int sf; // Switching function option.
#endif int nlevels_respa;
};
}
#endif
#endif

View File

@ -1,359 +1,365 @@
/* ------------------------------------------------------------------------- /* -------------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------- /* -------------------------------------------------------------------------
Contributing authors: Contributing authors:
Rodrigo Freitas (Unicamp/Brazil) - rodrigohb@gmail.com Rodrigo Freitas (Unicamp/Brazil) - rodrigohb@gmail.com
Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "stdlib.h" #include "stdlib.h"
#include "string.h" #include "string.h"
#include "fix_ti_spring.h" #include "fix_ti_spring.h"
#include "atom.h" #include "atom.h"
#include "update.h" #include "update.h"
#include "domain.h" #include "domain.h"
#include "respa.h" #include "respa.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
using namespace FixConst; using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
FixTISpring::FixTISpring(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) {
FixTISpring::FixTISpring(LAMMPS *lmp, int narg, char **arg) :
if (narg < 6 || narg > 8) Fix(lmp, narg, arg)
error->all(FLERR,"Illegal fix ti/spring command"); {
if (narg < 6 || narg > 8)
// Flags. error->all(FLERR,"Illegal fix ti/spring command");
restart_peratom = 1;
scalar_flag = 1; // Flags.
vector_flag = 1; restart_peratom = 1;
size_vector = 2; scalar_flag = 1;
global_freq = 1; global_freq = 1;
extscalar = 1; vector_flag = 1;
extvector = 1; size_vector = 2;
global_freq = 1;
// Spring constant. extscalar = 1;
k = atof(arg[3]); extvector = 1;
espring = 0.0;
if (k <= 0.0) error->all(FLERR,"Illegal fix ti/spring command"); // Spring constant.
k = force->numeric(FLERR,arg[3]);
// Initial position. if (k <= 0.0) error->all(FLERR,"Illegal fix ti/spring command");
xoriginal = NULL;
grow_arrays(atom->nmax); // Perform initial allocation of atom-based array.
atom->add_callback(0); // Registar with Atom class.
atom->add_callback(1); xoriginal = NULL;
grow_arrays(atom->nmax);
double **x = atom->x; atom->add_callback(0);
int *mask = atom->mask; atom->add_callback(1);
tagint *image = atom->image;
int nlocal = atom->nlocal; // xoriginal = initial unwrapped positions of atom.
for (int i = 0; i < nlocal; i++) { double **x = atom->x;
if (mask[i] & groupbit) domain->unmap(x[i], image[i], xoriginal[i]); int *mask = atom->mask;
else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; tagint *image = atom->image;
} int nlocal = atom->nlocal;
// Time variables. for (int i = 0; i < nlocal; i++) {
t_switch = atof(arg[4]); // Switching time. if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
t_equil = atof(arg[5]); // Equilibration time. else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
t0 = update->ntimestep; // Initial time. }
if (t_switch < 0.0) error->all(FLERR,"Illegal fix ti/spring command");
if (t_equil < 0.0) error->all(FLERR,"Illegal fix ti/spring command"); // Time variables.
t_switch = atoi(arg[4]); // Switching time.
// Coupling parameter initialization. t_equil = atoi(arg[5]); // Equilibration time.
sf = 1; t0 = update->ntimestep; // Initial time.
if (narg > 6) { if (t_switch < 0.0) error->all(FLERR,"Illegal fix ti/spring command");
if (strcmp(arg[6], "function") == 0) sf = atoi(arg[7]); if (t_equil < 0.0) error->all(FLERR,"Illegal fix ti/spring command");
else error->all(FLERR,"Illegal fix ti/spring switching function");
if ((sf!=1) && (sf!=2)) // Coupling parameter initialization.
error->all(FLERR,"Illegal fix ti/spring switching function"); sf = 1;
} if (narg > 6) {
lambda = switch_func(0); if (strcmp(arg[6], "function") == 0) sf = atoi(arg[7]);
dlambda = dswitch_func(0); else error->all(FLERR,"Illegal fix ti/spring switching function");
} if ((sf!=1) && (sf!=2))
error->all(FLERR,"Illegal fix ti/spring switching function");
/* ---------------------------------------------------------------------- */ }
lambda = switch_func(0);
FixTISpring::~FixTISpring() dlambda = dswitch_func(0);
{
// unregister callbacks to this fix from Atom class espring = 0.0;
atom->delete_callback(id,0); }
atom->delete_callback(id,1);
/* ---------------------------------------------------------------------- */
// delete locally stored array
memory->destroy(xoriginal); FixTISpring::~FixTISpring()
} {
// unregister callbacks to this fix from Atom class
/* ---------------------------------------------------------------------- */ atom->delete_callback(id,0);
atom->delete_callback(id,1);
int FixTISpring::setmask()
{ // delete locally stored array
int mask = 0; memory->destroy(xoriginal);
mask |= INITIAL_INTEGRATE; }
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA; /* ---------------------------------------------------------------------- */
mask |= MIN_POST_FORCE;
mask |= THERMO_ENERGY; int FixTISpring::setmask()
return mask; {
} int mask = 0;
mask |= INITIAL_INTEGRATE;
/* ---------------------------------------------------------------------- */ mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
void FixTISpring::init() mask |= MIN_POST_FORCE;
{ mask |= THERMO_ENERGY;
if (strstr(update->integrate_style,"respa")) return mask;
nlevels_respa = ((Respa *) update->integrate)->nlevels; }
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void FixTISpring::init()
void FixTISpring::setup(int vflag) {
{ if (strstr(update->integrate_style,"respa"))
if (strstr(update->integrate_style,"verlet")) nlevels_respa = ((Respa *) update->integrate)->nlevels;
post_force(vflag); }
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); /* ---------------------------------------------------------------------- */
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); void FixTISpring::setup(int vflag)
} {
} if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
/* ---------------------------------------------------------------------- */ else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
void FixTISpring::min_setup(int vflag) post_force_respa(vflag,nlevels_respa-1,0);
{ ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
post_force(vflag); }
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTISpring::post_force(int vflag) void FixTISpring::min_setup(int vflag)
{ {
post_force(vflag);
// If on the first equilibration do not calculate forces. }
int t = update->ntimestep - t0;
espring = 0.0; /* ---------------------------------------------------------------------- */
if(t < t_equil) return;
void FixTISpring::post_force(int vflag)
double **x = atom->x; {
double **f = atom->f; // If on the first equilibration do not calculate forces.
int *mask = atom->mask; int t = update->ntimestep - t0;
tagint *image = atom->image; if(t < t_equil) return;
int nlocal = atom->nlocal;
double **x = atom->x;
double dx, dy, dz; double **f = atom->f;
double unwrap[3]; int *mask = atom->mask;
tagint *image = atom->image;
for (int i = 0; i < nlocal; i++) int nlocal = atom->nlocal;
if (mask[i] & groupbit) {
domain->unmap(x[i], image[i], unwrap); double dx, dy, dz;
dx = unwrap[0] - xoriginal[i][0]; double unwrap[3];
dy = unwrap[1] - xoriginal[i][1];
dz = unwrap[2] - xoriginal[i][2]; espring = 0.0;
f[i][0] = (1-lambda) * f[i][0] + lambda * (-k*dx);
f[i][1] = (1-lambda) * f[i][1] + lambda * (-k*dy); for (int i = 0; i < nlocal; i++)
f[i][2] = (1-lambda) * f[i][2] + lambda * (-k*dz); if (mask[i] & groupbit) {
espring += k * (dx*dx + dy*dy + dz*dz); domain->unmap(x[i],image[i],unwrap);
} dx = unwrap[0] - xoriginal[i][0];
dy = unwrap[1] - xoriginal[i][1];
espring *= 0.5; dz = unwrap[2] - xoriginal[i][2];
} f[i][0] = (1-lambda) * f[i][0] + lambda * (-k*dx);
f[i][1] = (1-lambda) * f[i][1] + lambda * (-k*dy);
/* ---------------------------------------------------------------------- */ f[i][2] = (1-lambda) * f[i][2] + lambda * (-k*dz);
espring += k * (dx*dx + dy*dy + dz*dz);
void FixTISpring::post_force_respa(int vflag, int ilevel, int iloop) }
{
if (ilevel == nlevels_respa-1) post_force(vflag); espring *= 0.5;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTISpring::min_post_force(int vflag) void FixTISpring::post_force_respa(int vflag, int ilevel, int iloop)
{ {
post_force(vflag); if (ilevel == nlevels_respa-1) post_force(vflag);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTISpring::initial_integrate(int vflag) void FixTISpring::min_post_force(int vflag)
{ {
// Update the coupling parameter value. post_force(vflag);
double t = update->ntimestep - (t0+t_equil); }
if( (t >= 0) && (t <= t_switch) ) { /* ---------------------------------------------------------------------- */
lambda = switch_func(t/t_switch);
dlambda = dswitch_func(t/t_switch); void FixTISpring::initial_integrate(int vflag)
} {
// Update the coupling parameter value.
if( (t >= t_equil+t_switch) && (t <= (t_equil+2*t_switch)) ) { double t = update->ntimestep - (t0+t_equil);
lambda = switch_func(1.0 - (t - t_switch - t_equil)/t_switch );
dlambda = - dswitch_func(1.0 - (t - t_switch - t_equil)/t_switch ); if( (t >= 0) && (t <= t_switch) ) {
} lambda = switch_func(t/t_switch);
} dlambda = dswitch_func(t/t_switch);
}
/* ----------------------------------------------------------------------
energy of stretched springs if( (t >= t_equil+t_switch) && (t <= (t_equil+2*t_switch)) ) {
------------------------------------------------------------------------- */ lambda = switch_func(1.0 - (t - t_switch - t_equil)/t_switch );
dlambda = - dswitch_func(1.0 - (t - t_switch - t_equil)/t_switch );
double FixTISpring::compute_scalar() }
{ }
double all;
MPI_Allreduce(&espring,&all,1,MPI_DOUBLE,MPI_SUM,world); /* ----------------------------------------------------------------------
return all; energy of stretched springs
} ------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- double FixTISpring::compute_scalar()
information about coupling parameter {
------------------------------------------------------------------------- */ double all;
MPI_Allreduce(&espring,&all,1,MPI_DOUBLE,MPI_SUM,world);
double FixTISpring::compute_vector(int n) return all;
{ }
linfo[0] = lambda;
linfo[1] = dlambda; /* ----------------------------------------------------------------------
return linfo[n]; information about coupling parameter
} ------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- double FixTISpring::compute_vector(int n)
memory usage of local atom-based array {
------------------------------------------------------------------------- */ linfo[0] = lambda;
linfo[1] = dlambda;
double FixTISpring::memory_usage() return linfo[n];
{ }
double bytes = atom->nmax*3 * sizeof(double);
return bytes; /* ----------------------------------------------------------------------
} memory usage of local atom-based array
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
allocate atom-based array double FixTISpring::memory_usage()
------------------------------------------------------------------------- */ {
double bytes = atom->nmax*3 * sizeof(double);
void FixTISpring::grow_arrays(int nmax) return bytes;
{ }
memory->grow(xoriginal,nmax,3,"fix_ti/spring:xoriginal");
} /* ----------------------------------------------------------------------
allocate atom-based array
/* ---------------------------------------------------------------------- ------------------------------------------------------------------------- */
copy values within local atom-based array
------------------------------------------------------------------------- */ void FixTISpring::grow_arrays(int nmax)
{
void FixTISpring::copy_arrays(int i, int j) memory->grow(xoriginal,nmax,3,"fix_ti/spring:xoriginal");
{ }
xoriginal[j][0] = xoriginal[i][0];
xoriginal[j][1] = xoriginal[i][1]; /* ----------------------------------------------------------------------
xoriginal[j][2] = xoriginal[i][2]; copy values within local atom-based array
} ------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- void FixTISpring::copy_arrays(int i, int j, int delflag)
pack values in local atom-based array for exchange with another proc {
------------------------------------------------------------------------- */ xoriginal[j][0] = xoriginal[i][0];
xoriginal[j][1] = xoriginal[i][1];
int FixTISpring::pack_exchange(int i, double *buf) xoriginal[j][2] = xoriginal[i][2];
{ }
buf[0] = xoriginal[i][0];
buf[1] = xoriginal[i][1]; /* ----------------------------------------------------------------------
buf[2] = xoriginal[i][2]; pack values in local atom-based array for exchange with another proc
return 3; ------------------------------------------------------------------------- */
}
int FixTISpring::pack_exchange(int i, double *buf)
/* ---------------------------------------------------------------------- {
unpack values in local atom-based array from exchange with another proc buf[0] = xoriginal[i][0];
------------------------------------------------------------------------- */ buf[1] = xoriginal[i][1];
buf[2] = xoriginal[i][2];
int FixTISpring::unpack_exchange(int nlocal, double *buf) return 3;
{ }
xoriginal[nlocal][0] = buf[0];
xoriginal[nlocal][1] = buf[1]; /* ----------------------------------------------------------------------
xoriginal[nlocal][2] = buf[2]; unpack values in local atom-based array from exchange with another proc
return 3; ------------------------------------------------------------------------- */
}
int FixTISpring::unpack_exchange(int nlocal, double *buf)
/* ---------------------------------------------------------------------- {
pack values in local atom-based arrays for restart file xoriginal[nlocal][0] = buf[0];
------------------------------------------------------------------------- */ xoriginal[nlocal][1] = buf[1];
xoriginal[nlocal][2] = buf[2];
int FixTISpring::pack_restart(int i, double *buf) return 3;
{ }
buf[0] = 4;
buf[1] = xoriginal[i][0]; /* ----------------------------------------------------------------------
buf[2] = xoriginal[i][1]; pack values in local atom-based arrays for restart file
buf[3] = xoriginal[i][2]; ------------------------------------------------------------------------- */
return 4;
} int FixTISpring::pack_restart(int i, double *buf)
{
buf[0] = 4;
/* ---------------------------------------------------------------------- buf[1] = xoriginal[i][0];
unpack values from atom->extra array to restart the fix buf[2] = xoriginal[i][1];
------------------------------------------------------------------------- */ buf[3] = xoriginal[i][2];
return 4;
void FixTISpring::unpack_restart(int nlocal, int nth) }
{
double **extra = atom->extra; /* ----------------------------------------------------------------------
unpack values from atom->extra array to restart the fix
// skip to Nth set of extra values ------------------------------------------------------------------------- */
int m = 0; void FixTISpring::unpack_restart(int nlocal, int nth)
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]); {
m++; double **extra = atom->extra;
xoriginal[nlocal][0] = extra[nlocal][m++]; // skip to Nth set of extra values
xoriginal[nlocal][1] = extra[nlocal][m++];
xoriginal[nlocal][2] = extra[nlocal][m++]; int m = 0;
} for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;
/* ----------------------------------------------------------------------
maxsize of any atom's restart data xoriginal[nlocal][0] = extra[nlocal][m++];
------------------------------------------------------------------------- */ xoriginal[nlocal][1] = extra[nlocal][m++];
xoriginal[nlocal][2] = extra[nlocal][m++];
int FixTISpring::maxsize_restart() }
{
return 4; /* ----------------------------------------------------------------------
} maxsize of any atom's restart data
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
size of atom nlocal's restart data int FixTISpring::maxsize_restart()
------------------------------------------------------------------------- */ {
return 4;
int FixTISpring::size_restart(int nlocal) }
{
return 4; /* ----------------------------------------------------------------------
} size of atom nlocal's restart data
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Switching function. int FixTISpring::size_restart(int nlocal)
------------------------------------------------------------------------- */ {
return 4;
double FixTISpring::switch_func(double t) }
{
if (sf == 1) return t; /* ----------------------------------------------------------------------
Switching function.
double t2 = t*t; ------------------------------------------------------------------------- */
double t5 = t2*t2*t;
return ((70.0*t2*t2 - 315.0*t2*t + 540.0*t2 - 420.0*t + 126.0)*t5); double FixTISpring::switch_func(double t)
} {
if (sf == 1) return t;
/* ----------------------------------------------------------------------
Switching function derivative. double t2 = t*t;
------------------------------------------------------------------------- */ double t5 = t2*t2*t;
return ((70.0*t2*t2 - 315.0*t2*t + 540.0*t2 - 420.0*t + 126.0)*t5);
double FixTISpring::dswitch_func(double t) }
{
if(sf == 1) return 1.0/t_switch; /* ----------------------------------------------------------------------
Switching function derivative.
double t2 = t*t; ------------------------------------------------------------------------- */
double t4 = t2*t2;
return ((630*t2*t2 - 2520*t2*t + 3780*t2 - 2520*t + 630)*t4) / t_switch; double FixTISpring::dswitch_func(double t)
} {
if(sf == 1) return 1.0/t_switch;
double t2 = t*t;
double t4 = t2*t2;
return ((630*t2*t2 - 2520*t2*t + 3780*t2 - 2520*t + 630)*t4) / t_switch;
}

View File

@ -1,88 +1,94 @@
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#ifdef FIX_CLASS /* ----------------------------------------------------------------------
Contributing authors:
FixStyle(ti/spring,FixTISpring) Rodrigo Freitas (Unicamp/Brazil) - rodrigohb@gmail.com
Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br
#else ------------------------------------------------------------------------- */
#ifndef LMP_FIX_TI_SPRING_H #ifdef FIX_CLASS
#define LMP_FIX_TI_SPRING_H
FixStyle(ti/spring,FixTISpring)
#include "fix.h"
#else
namespace LAMMPS_NS {
#ifndef LMP_FIX_TI_SPRING_H
class FixTISpring : public Fix { #define LMP_FIX_TI_SPRING_H
public:
FixTISpring(class LAMMPS *, int, char **); #include "fix.h"
~FixTISpring();
int setmask(); namespace LAMMPS_NS {
void init();
void setup(int); class FixTISpring : public Fix {
void min_setup(int); public:
void post_force(int); FixTISpring(class LAMMPS *, int, char **);
void post_force_respa(int, int, int); ~FixTISpring();
void min_post_force(int); int setmask();
void initial_integrate(int); void init();
double compute_scalar(); void setup(int);
double compute_vector(int); void min_setup(int);
void post_force(int);
double memory_usage(); void post_force_respa(int, int, int);
void grow_arrays(int); void min_post_force(int);
void copy_arrays(int, int); void initial_integrate(int);
int pack_exchange(int, double *); double compute_scalar();
int unpack_exchange(int, double *); double compute_vector(int);
int pack_restart(int, double *);
void unpack_restart(int, int); double memory_usage();
int size_restart(int); void grow_arrays(int);
int maxsize_restart(); void copy_arrays(int, int, int);
int pack_exchange(int, double *);
private: int unpack_exchange(int, double *);
double switch_func(double); // Switching function. int pack_restart(int, double *);
double dswitch_func(double); // Switching function derivative. void unpack_restart(int, int);
int size_restart(int);
double k; // Spring constant. int maxsize_restart();
double espring; // Springs energies.
double **xoriginal; // Original coords of atoms. private:
double lambda; // Coupling parameter. double switch_func(double); // Switching function.
double dlambda; // Lambda variation with t. double dswitch_func(double); // Switching function derivative.
double linfo[2]; // Current lambda status.
int t_switch; // Total switching steps. double k; // Spring constant.
int t_equil; // Equilibration time. double espring; // Springs energies.
int t0; // Initial time. double **xoriginal; // Original coords of atoms.
int sf; // Switching function option. double lambda; // Coupling parameter.
int nlevels_respa; double dlambda; // Lambda variation with t.
}; double linfo[2]; // Current lambda status.
int t_switch; // Total switching steps.
} int t_equil; // Equilibration time.
int t0; // Initial time.
#endif int sf; // Switching function option.
#endif int nlevels_respa;
};
/* ERROR/WARNING messages:
}
E: Illegal ... command
#endif
Self-explanatory. Check the input script syntax and compare to the #endif
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line. /* ERROR/WARNING messages:
E: Illegal fix ti/spring switching function E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line. command-line option when running LAMMPS to see the offending line.
*/ E: Illegal fix ti/spring switching function
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.
*/