ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2013-11-07 14:38:21 +00:00
parent 6c76b65621
commit 05649f1189
4 changed files with 722 additions and 722 deletions
Download Patch File Download Diff File Expand all files Collapse all files

394
src/USER-MISC/fix_ti_rs.cpp
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 = atoi(arg[5]); t_switch = atoi(arg[5]);
t_equil = atoi(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;
} }

132
src/USER-MISC/fix_ti_rs.h
View File

@ -1,66 +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.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
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
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#ifdef FIX_CLASS #ifdef FIX_CLASS
FixStyle(ti/rs,FixTIRS) FixStyle(ti/rs,FixTIRS)
#else #else
#ifndef LMP_FIX_TI_RS_H #ifndef LMP_FIX_TI_RS_H
#define LMP_FIX_TI_RS_H #define LMP_FIX_TI_RS_H
#include "fix.h" #include "fix.h"
namespace LAMMPS_NS { namespace LAMMPS_NS {
class FixTIRS : public Fix { class FixTIRS : public Fix {
public: public:
FixTIRS(class LAMMPS *, int, char **); FixTIRS(class LAMMPS *, int, char **);
~FixTIRS(); ~FixTIRS();
int setmask(); int setmask();
void init(); void init();
void setup(int); void setup(int);
void min_setup(int); void min_setup(int);
void post_force(int); void post_force(int);
void post_force_respa(int, int, int); void post_force_respa(int, int, int);
void min_post_force(int); void min_post_force(int);
virtual void initial_integrate(int); virtual void initial_integrate(int);
double compute_vector(int); double compute_vector(int);
private: private:
double switch_func(double); double switch_func(double);
double dswitch_func(double); double dswitch_func(double);
double lambda; // Coupling parameter. double lambda; // Coupling parameter.
double dlambda; // Lambda variation with t. double dlambda; // Lambda variation with t.
double l_initial; // Lambda initial value. double l_initial; // Lambda initial value.
double l_final; // Lambda final value. double l_final; // Lambda final value.
double linfo[2]; // Current lambda status. double linfo[2]; // Current lambda status.
int t_switch; // Total switching steps. int t_switch; // Total switching steps.
int t_equil; // Equilibration time. int t_equil; // Equilibration time.
int t0; // Initial time. int t0; // Initial time.
int sf; // Switching function option. int sf; // Switching function option.
int nlevels_respa; int nlevels_respa;
}; };
} }
#endif #endif
#endif #endif

730
src/USER-MISC/fix_ti_spring.cpp
View File

@ -1,365 +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" #include "force.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixTISpring::FixTISpring(LAMMPS *lmp, int narg, char **arg) : FixTISpring::FixTISpring(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) Fix(lmp, narg, arg)
{ {
if (narg < 6 || narg > 8) if (narg < 6 || narg > 8)
error->all(FLERR,"Illegal fix ti/spring command"); error->all(FLERR,"Illegal fix ti/spring command");
// Flags. // Flags.
restart_peratom = 1; restart_peratom = 1;
scalar_flag = 1; scalar_flag = 1;
global_freq = 1; global_freq = 1;
vector_flag = 1; vector_flag = 1;
size_vector = 2; size_vector = 2;
global_freq = 1; global_freq = 1;
extscalar = 1; extscalar = 1;
extvector = 1; extvector = 1;
// Spring constant. // Spring constant.
k = force->numeric(FLERR,arg[3]); k = force->numeric(FLERR,arg[3]);
if (k <= 0.0) error->all(FLERR,"Illegal fix ti/spring command"); if (k <= 0.0) error->all(FLERR,"Illegal fix ti/spring command");
// Perform initial allocation of atom-based array. // Perform initial allocation of atom-based array.
// Registar with Atom class. // Registar with Atom class.
xoriginal = NULL; xoriginal = NULL;
grow_arrays(atom->nmax); grow_arrays(atom->nmax);
atom->add_callback(0); atom->add_callback(0);
atom->add_callback(1); atom->add_callback(1);
// xoriginal = initial unwrapped positions of atom. // xoriginal = initial unwrapped positions of atom.
double **x = atom->x; double **x = atom->x;
int *mask = atom->mask; int *mask = atom->mask;
tagint *image = atom->image; tagint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]); if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
} }
// Time variables. // Time variables.
t_switch = atoi(arg[4]); // Switching time. t_switch = atoi(arg[4]); // Switching time.
t_equil = atoi(arg[5]); // Equilibration time. t_equil = atoi(arg[5]); // Equilibration time.
t0 = update->ntimestep; // Initial time. t0 = update->ntimestep; // Initial time.
if (t_switch < 0.0) error->all(FLERR,"Illegal fix ti/spring command"); 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"); if (t_equil < 0.0) error->all(FLERR,"Illegal fix ti/spring command");
// Coupling parameter initialization. // Coupling parameter initialization.
sf = 1; sf = 1;
if (narg > 6) { if (narg > 6) {
if (strcmp(arg[6], "function") == 0) sf = atoi(arg[7]); if (strcmp(arg[6], "function") == 0) sf = atoi(arg[7]);
else error->all(FLERR,"Illegal fix ti/spring switching function"); else error->all(FLERR,"Illegal fix ti/spring switching function");
if ((sf!=1) && (sf!=2)) if ((sf!=1) && (sf!=2))
error->all(FLERR,"Illegal fix ti/spring switching function"); error->all(FLERR,"Illegal fix ti/spring switching function");
} }
lambda = switch_func(0); lambda = switch_func(0);
dlambda = dswitch_func(0); dlambda = dswitch_func(0);
espring = 0.0; espring = 0.0;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixTISpring::~FixTISpring() FixTISpring::~FixTISpring()
{ {
// 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);
// delete locally stored array // delete locally stored array
memory->destroy(xoriginal); memory->destroy(xoriginal);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
int FixTISpring::setmask() int FixTISpring::setmask()
{ {
int mask = 0; int mask = 0;
mask |= INITIAL_INTEGRATE; mask |= INITIAL_INTEGRATE;
mask |= POST_FORCE; mask |= POST_FORCE;
mask |= POST_FORCE_RESPA; mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE; mask |= MIN_POST_FORCE;
mask |= THERMO_ENERGY; mask |= THERMO_ENERGY;
return mask; return mask;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTISpring::init() void FixTISpring::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 FixTISpring::setup(int vflag) void FixTISpring::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 FixTISpring::min_setup(int vflag) void FixTISpring::min_setup(int vflag)
{ {
post_force(vflag); post_force(vflag);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTISpring::post_force(int vflag) void FixTISpring::post_force(int vflag)
{ {
// If on the first equilibration do not calculate forces. // If on the first equilibration do not calculate forces.
int t = update->ntimestep - t0; int t = update->ntimestep - t0;
if(t < t_equil) return; if(t < t_equil) return;
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
int *mask = atom->mask; int *mask = atom->mask;
tagint *image = atom->image; tagint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double dx, dy, dz; double dx, dy, dz;
double unwrap[3]; double unwrap[3];
espring = 0.0; espring = 0.0;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xoriginal[i][0]; dx = unwrap[0] - xoriginal[i][0];
dy = unwrap[1] - xoriginal[i][1]; dy = unwrap[1] - xoriginal[i][1];
dz = unwrap[2] - xoriginal[i][2]; dz = unwrap[2] - xoriginal[i][2];
f[i][0] = (1-lambda) * f[i][0] + lambda * (-k*dx); f[i][0] = (1-lambda) * f[i][0] + lambda * (-k*dx);
f[i][1] = (1-lambda) * f[i][1] + lambda * (-k*dy); f[i][1] = (1-lambda) * f[i][1] + lambda * (-k*dy);
f[i][2] = (1-lambda) * f[i][2] + lambda * (-k*dz); f[i][2] = (1-lambda) * f[i][2] + lambda * (-k*dz);
espring += k * (dx*dx + dy*dy + dz*dz); espring += k * (dx*dx + dy*dy + dz*dz);
} }
espring *= 0.5; espring *= 0.5;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTISpring::post_force_respa(int vflag, int ilevel, int iloop) void FixTISpring::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 FixTISpring::min_post_force(int vflag) void FixTISpring::min_post_force(int vflag)
{ {
post_force(vflag); post_force(vflag);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTISpring::initial_integrate(int vflag) void FixTISpring::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 );
} }
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
energy of stretched springs energy of stretched springs
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
double FixTISpring::compute_scalar() double FixTISpring::compute_scalar()
{ {
double all; double all;
MPI_Allreduce(&espring,&all,1,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(&espring,&all,1,MPI_DOUBLE,MPI_SUM,world);
return all; return all;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
information about coupling parameter information about coupling parameter
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
double FixTISpring::compute_vector(int n) double FixTISpring::compute_vector(int n)
{ {
linfo[0] = lambda; linfo[0] = lambda;
linfo[1] = dlambda; linfo[1] = dlambda;
return linfo[n]; return linfo[n];
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
memory usage of local atom-based array memory usage of local atom-based array
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
double FixTISpring::memory_usage() double FixTISpring::memory_usage()
{ {
double bytes = atom->nmax*3 * sizeof(double); double bytes = atom->nmax*3 * sizeof(double);
return bytes; return bytes;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
allocate atom-based array allocate atom-based array
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixTISpring::grow_arrays(int nmax) void FixTISpring::grow_arrays(int nmax)
{ {
memory->grow(xoriginal,nmax,3,"fix_ti/spring:xoriginal"); memory->grow(xoriginal,nmax,3,"fix_ti/spring:xoriginal");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
copy values within local atom-based array copy values within local atom-based array
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixTISpring::copy_arrays(int i, int j, int delflag) void FixTISpring::copy_arrays(int i, int j, int delflag)
{ {
xoriginal[j][0] = xoriginal[i][0]; xoriginal[j][0] = xoriginal[i][0];
xoriginal[j][1] = xoriginal[i][1]; xoriginal[j][1] = xoriginal[i][1];
xoriginal[j][2] = xoriginal[i][2]; xoriginal[j][2] = xoriginal[i][2];
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int FixTISpring::pack_exchange(int i, double *buf) int FixTISpring::pack_exchange(int i, double *buf)
{ {
buf[0] = xoriginal[i][0]; buf[0] = xoriginal[i][0];
buf[1] = xoriginal[i][1]; buf[1] = xoriginal[i][1];
buf[2] = xoriginal[i][2]; buf[2] = xoriginal[i][2];
return 3; return 3;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int FixTISpring::unpack_exchange(int nlocal, double *buf) int FixTISpring::unpack_exchange(int nlocal, double *buf)
{ {
xoriginal[nlocal][0] = buf[0]; xoriginal[nlocal][0] = buf[0];
xoriginal[nlocal][1] = buf[1]; xoriginal[nlocal][1] = buf[1];
xoriginal[nlocal][2] = buf[2]; xoriginal[nlocal][2] = buf[2];
return 3; return 3;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
pack values in local atom-based arrays for restart file pack values in local atom-based arrays for restart file
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int FixTISpring::pack_restart(int i, double *buf) int FixTISpring::pack_restart(int i, double *buf)
{ {
buf[0] = 4; buf[0] = 4;
buf[1] = xoriginal[i][0]; buf[1] = xoriginal[i][0];
buf[2] = xoriginal[i][1]; buf[2] = xoriginal[i][1];
buf[3] = xoriginal[i][2]; buf[3] = xoriginal[i][2];
return 4; return 4;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
unpack values from atom->extra array to restart the fix unpack values from atom->extra array to restart the fix
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixTISpring::unpack_restart(int nlocal, int nth) void FixTISpring::unpack_restart(int nlocal, int nth)
{ {
double **extra = atom->extra; double **extra = atom->extra;
// skip to Nth set of extra values // skip to Nth set of extra values
int m = 0; int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]); for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++; m++;
xoriginal[nlocal][0] = extra[nlocal][m++]; xoriginal[nlocal][0] = extra[nlocal][m++];
xoriginal[nlocal][1] = extra[nlocal][m++]; xoriginal[nlocal][1] = extra[nlocal][m++];
xoriginal[nlocal][2] = extra[nlocal][m++]; xoriginal[nlocal][2] = extra[nlocal][m++];
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
maxsize of any atom's restart data maxsize of any atom's restart data
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int FixTISpring::maxsize_restart() int FixTISpring::maxsize_restart()
{ {
return 4; return 4;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
size of atom nlocal's restart data size of atom nlocal's restart data
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int FixTISpring::size_restart(int nlocal) int FixTISpring::size_restart(int nlocal)
{ {
return 4; return 4;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
Switching function. Switching function.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
double FixTISpring::switch_func(double t) double FixTISpring::switch_func(double t)
{ {
if (sf == 1) return t; if (sf == 1) return t;
double t2 = t*t; double t2 = t*t;
double t5 = t2*t2*t; double t5 = t2*t2*t;
return ((70.0*t2*t2 - 315.0*t2*t + 540.0*t2 - 420.0*t + 126.0)*t5); return ((70.0*t2*t2 - 315.0*t2*t + 540.0*t2 - 420.0*t + 126.0)*t5);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
Switching function derivative. Switching function derivative.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
double FixTISpring::dswitch_func(double t) double FixTISpring::dswitch_func(double t)
{ {
if(sf == 1) return 1.0/t_switch; if(sf == 1) return 1.0/t_switch;
double t2 = t*t; double t2 = t*t;
double t4 = t2*t2; double t4 = t2*t2;
return ((630*t2*t2 - 2520*t2*t + 3780*t2 - 2520*t + 630)*t4) / t_switch; return ((630*t2*t2 - 2520*t2*t + 3780*t2 - 2520*t + 630)*t4) / t_switch;
} }

188
src/USER-MISC/fix_ti_spring.h
View File

@ -1,94 +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.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
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
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#ifdef FIX_CLASS #ifdef FIX_CLASS
FixStyle(ti/spring,FixTISpring) FixStyle(ti/spring,FixTISpring)
#else #else
#ifndef LMP_FIX_TI_SPRING_H #ifndef LMP_FIX_TI_SPRING_H
#define LMP_FIX_TI_SPRING_H #define LMP_FIX_TI_SPRING_H
#include "fix.h" #include "fix.h"
namespace LAMMPS_NS { namespace LAMMPS_NS {
class FixTISpring : public Fix { class FixTISpring : public Fix {
public: public:
FixTISpring(class LAMMPS *, int, char **); FixTISpring(class LAMMPS *, int, char **);
~FixTISpring(); ~FixTISpring();
int setmask(); int setmask();
void init(); void init();
void setup(int); void setup(int);
void min_setup(int); void min_setup(int);
void post_force(int); void post_force(int);
void post_force_respa(int, int, int); void post_force_respa(int, int, int);
void min_post_force(int); void min_post_force(int);
void initial_integrate(int); void initial_integrate(int);
double compute_scalar(); double compute_scalar();
double compute_vector(int); double compute_vector(int);
double memory_usage(); double memory_usage();
void grow_arrays(int); void grow_arrays(int);
void copy_arrays(int, int, int); void copy_arrays(int, int, int);
int pack_exchange(int, double *); int pack_exchange(int, double *);
int unpack_exchange(int, double *); int unpack_exchange(int, double *);
int pack_restart(int, double *); int pack_restart(int, double *);
void unpack_restart(int, int); void unpack_restart(int, int);
int size_restart(int); int size_restart(int);
int maxsize_restart(); int maxsize_restart();
private: private:
double switch_func(double); // Switching function. double switch_func(double); // Switching function.
double dswitch_func(double); // Switching function derivative. double dswitch_func(double); // Switching function derivative.
double k; // Spring constant. double k; // Spring constant.
double espring; // Springs energies. double espring; // Springs energies.
double **xoriginal; // Original coords of atoms. double **xoriginal; // Original coords of atoms.
double lambda; // Coupling parameter. double lambda; // Coupling parameter.
double dlambda; // Lambda variation with t. double dlambda; // Lambda variation with t.
double linfo[2]; // Current lambda status. double linfo[2]; // Current lambda status.
int t_switch; // Total switching steps. int t_switch; // Total switching steps.
int t_equil; // Equilibration time. int t_equil; // Equilibration time.
int t0; // Initial time. int t0; // Initial time.
int sf; // Switching function option. int sf; // Switching function option.
int nlevels_respa; int nlevels_respa;
}; };
} }
#endif #endif
#endif #endif
/* ERROR/WARNING messages: /* ERROR/WARNING messages:
E: Illegal ... command 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 E: Illegal fix ti/spring switching function
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.
*/ */