diff --git a/src/USER-MISC/fix_ti_rs.cpp b/src/USER-MISC/fix_ti_rs.cpp deleted file mode 100755 index 1090fc02d..000000000 --- a/src/USER-MISC/fix_ti_rs.cpp +++ /dev/null @@ -1,199 +0,0 @@ -/* ------------------------------------------------------------------------- - 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 authors: - Rodrigo Freitas (Unicamp/Brazil) - rodrigohb@gmail.com - Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br -------------------------------------------------------------------------- */ - -#include <stdlib.h> -#include <string.h> -#include <math.h> -#include "fix_ti_rs.h" -#include "atom.h" -#include "update.h" -#include "respa.h" -#include "error.h" -#include "force.h" - -using namespace LAMMPS_NS; -using namespace FixConst; - -/* ---------------------------------------------------------------------- */ - -// Class constructor initialize all variables. - -FixTIRS::FixTIRS(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) -{ - // Checking the input information. - if (narg < 7 || narg > 9) error->all(FLERR,"Illegal fix ti/rs command"); - - // Fix flags. - vector_flag = 1; - size_vector = 2; - global_freq = 1; - extvector = 1; - - // Time variables. - t_switch = force->bnumeric(FLERR,arg[5]); - t_equil = force->bnumeric(FLERR,arg[6]); - t0 = update->ntimestep; - if (t_switch <= 0) error->all(FLERR,"Illegal fix ti/rs command"); - if (t_equil <= 0) error->all(FLERR,"Illegal fix ti/rs command"); - - // Coupling parameter limits and initialization. - l_initial = force->numeric(FLERR,arg[3]); - l_final = force->numeric(FLERR,arg[4]); - sf = 1; - if (narg > 7) { - if (strcmp(arg[7], "function") == 0) sf = force->inumeric(FLERR,arg[8]); - else error->all(FLERR,"Illegal fix ti/rs switching function"); - if ((sf<1) || (sf>3)) - error->all(FLERR,"Illegal fix ti/rs switching function"); - } - lambda = switch_func(0); - dlambda = dswitch_func(0); -} - -/* ---------------------------------------------------------------------- */ - -FixTIRS::~FixTIRS() -{ - // unregister callbacks to this fix from Atom class - atom->delete_callback(id,0); - atom->delete_callback(id,1); -} - -/* ---------------------------------------------------------------------- */ - -int FixTIRS::setmask() -{ - int mask = 0; - mask |= INITIAL_INTEGRATE; - mask |= POST_FORCE; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixTIRS::init() -{ - if (strstr(update->integrate_style,"respa")) - nlevels_respa = ((Respa *) update->integrate)->nlevels; -} - -/* ---------------------------------------------------------------------- */ - -void FixTIRS::setup(int vflag) -{ - if (strstr(update->integrate_style,"verlet")) - 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 FixTIRS::min_setup(int vflag) -{ - post_force(vflag); -} - - -/* ---------------------------------------------------------------------- */ - -void FixTIRS::post_force(int vflag) -{ - - int *mask = atom->mask; - int nlocal = atom->nlocal; - double **f = atom->f; - - // Scaling forces. - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - f[i][0] = lambda * f[i][0]; - f[i][1] = lambda * f[i][1]; - f[i][2] = lambda * f[i][2]; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixTIRS::post_force_respa(int vflag, int ilevel, int iloop) -{ - if (ilevel == nlevels_respa-1) post_force(vflag); -} - -/* ---------------------------------------------------------------------- */ - -void FixTIRS::min_post_force(int vflag) -{ - post_force(vflag); -} - -/* ---------------------------------------------------------------------- */ - -void FixTIRS::initial_integrate(int vflag) -{ - // Update the coupling parameter value. - const bigint t = update->ntimestep - (t0+t_equil); - const double r_switch = 1.0/t_switch; - - if( (t >= 0) && (t <= t_switch) ) { - lambda = switch_func(t*r_switch); - dlambda = dswitch_func(t*r_switch); - } - - if( (t >= t_equil+t_switch) && (t <= (t_equil+2*t_switch)) ) { - lambda = switch_func(1.0 - (t - t_switch - t_equil)*r_switch); - dlambda = - dswitch_func(1.0 - (t - t_switch - t_equil)*r_switch); - } -} - -/* ---------------------------------------------------------------------- */ - -double FixTIRS::compute_vector(int n) -{ - linfo[0] = lambda; - linfo[1] = dlambda; - return linfo[n]; -} - -/* ---------------------------------------------------------------------- */ - -double FixTIRS::switch_func(double t) -{ - 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)); - - // Default option is sf = 1. - return l_initial + (l_final - l_initial) * t; -} - -/* ---------------------------------------------------------------------- */ - -double FixTIRS::dswitch_func(double t) -{ - double aux = (1.0/l_initial - 1.0/l_final); - if (sf == 2) return lambda * lambda * aux / t_switch; - if (sf == 3) return lambda * lambda * aux / (t_switch * log(2) * (1 + t)); - - // Default option is sf = 1. - return (l_final-l_initial)/t_switch; -} diff --git a/src/USER-MISC/fix_ti_rs.h b/src/USER-MISC/fix_ti_rs.h deleted file mode 100755 index 0aee21a87..000000000 --- a/src/USER-MISC/fix_ti_rs.h +++ /dev/null @@ -1,66 +0,0 @@ -/* -*- 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. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: - Rodrigo Freitas (Unicamp/Brazil) - rodrigohb@gmail.com - Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br -------------------------------------------------------------------------- */ - -#ifdef FIX_CLASS - -FixStyle(ti/rs,FixTIRS) - -#else - -#ifndef LMP_FIX_TI_RS_H -#define LMP_FIX_TI_RS_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixTIRS : public Fix { - public: - FixTIRS(class LAMMPS *, int, char **); - ~FixTIRS(); - int setmask(); - void init(); - void setup(int); - void min_setup(int); - void post_force(int); - void post_force_respa(int, int, int); - void min_post_force(int); - virtual void initial_integrate(int); - double compute_vector(int); - - private: - double switch_func(double); - double dswitch_func(double); - - double lambda; // Coupling parameter. - double dlambda; // Lambda variation with t. - double l_initial; // Lambda initial value. - double l_final; // Lambda final value. - double linfo[2]; // Current lambda status. - bigint t_switch; // Total switching steps. - bigint t_equil; // Equilibration time. - bigint t0; // Initial time. - int sf; // Switching function option. - int nlevels_respa; -}; - -} - -#endif -#endif diff --git a/src/USER-MISC/fix_ti_spring.cpp b/src/USER-MISC/fix_ti_spring.cpp index edb6a7b2c..5f3e9d9cd 100755 --- a/src/USER-MISC/fix_ti_spring.cpp +++ b/src/USER-MISC/fix_ti_spring.cpp @@ -1,366 +1,378 @@ /* ------------------------------------------------------------------------- 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 authors: - Rodrigo Freitas (Unicamp/Brazil) - rodrigohb@gmail.com + Contributing authors: + Rodrigo Freitas (UC Berkeley) - rodrigof@berkeley.edu + Mark Asta (UC Berkeley) - mdasta@berkeley.edu Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br ------------------------------------------------------------------------- */ -#include <stdlib.h> -#include <string.h> +#include "stdlib.h" +#include "string.h" #include "fix_ti_spring.h" #include "atom.h" #include "update.h" #include "domain.h" #include "respa.h" #include "memory.h" #include "error.h" #include "force.h" -using namespace LAMMPS_NS; -using namespace FixConst; +using namespace LAMMPS_NS; +using namespace FixConst; + +static const char cite_ti_spring[] = + "ti/spring command:\n\n" + "@article{freitas2016,\n" + " author={Freitas, Rodrigo and Asta, Mark and de Koning, Maurice},\n" + " title={Nonequilibrium free-energy calculation of solids using LAMMPS},\n" + " journal={Computational Materials Science},\n" + " volume={112},\n" + " pages={333--341},\n" + " year={2016},\n" + " publisher={Elsevier}\n" + "}\n\n"; /* ---------------------------------------------------------------------- */ -FixTISpring::FixTISpring(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) +FixTISpring::FixTISpring(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) { if (narg < 6 || narg > 8) error->all(FLERR,"Illegal fix ti/spring command"); // Flags. restart_peratom = 1; scalar_flag = 1; global_freq = 1; vector_flag = 1; size_vector = 2; global_freq = 1; extscalar = 1; extvector = 1; // Spring constant. k = force->numeric(FLERR,arg[3]); if (k <= 0.0) error->all(FLERR,"Illegal fix ti/spring command"); // Perform initial allocation of atom-based array. // Registar with Atom class. xoriginal = NULL; grow_arrays(atom->nmax); atom->add_callback(0); atom->add_callback(1); // xoriginal = initial unwrapped positions of atom. double **x = atom->x; int *mask = atom->mask; - imageint *image = atom->image; + tagint *image = atom->image; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; 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; } // Time variables. - t_switch = force->bnumeric(FLERR,arg[4]); // Number of steps for switching. - t_equil = force->bnumeric(FLERR,arg[5]); // Number of steps for equilibration. + t_switch = atoi(arg[4]); // Switching time. + t_equil = atoi(arg[5]); // Equilibration time. t0 = update->ntimestep; // Initial time. - if (t_switch <= 0) error->all(FLERR,"Illegal fix ti/spring command"); - if (t_equil <= 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"); // Coupling parameter initialization. sf = 1; if (narg > 6) { - if (strcmp(arg[6], "function") == 0) sf = force->inumeric(FLERR,arg[7]); + if (strcmp(arg[6], "function") == 0) sf = atoi(arg[7]); 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"); } - lambda = switch_func(0); - dlambda = dswitch_func(0); + lambda = switch_func(0); + dlambda = dswitch_func(0); espring = 0.0; } /* ---------------------------------------------------------------------- */ FixTISpring::~FixTISpring() { // unregister callbacks to this fix from Atom class atom->delete_callback(id,0); atom->delete_callback(id,1); // delete locally stored array memory->destroy(xoriginal); } /* ---------------------------------------------------------------------- */ int FixTISpring::setmask() { int mask = 0; mask |= INITIAL_INTEGRATE; mask |= POST_FORCE; mask |= POST_FORCE_RESPA; mask |= MIN_POST_FORCE; mask |= THERMO_ENERGY; return mask; } /* ---------------------------------------------------------------------- */ void FixTISpring::init() { if (strstr(update->integrate_style,"respa")) nlevels_respa = ((Respa *) update->integrate)->nlevels; } /* ---------------------------------------------------------------------- */ void FixTISpring::setup(int vflag) { if (strstr(update->integrate_style,"verlet")) 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::min_setup(int vflag) { post_force(vflag); } /* ---------------------------------------------------------------------- */ void FixTISpring::post_force(int vflag) { // If on the first equilibration do not calculate forces. - bigint t = update->ntimestep - t0; + int t = update->ntimestep - t0; if(t < t_equil) return; double **x = atom->x; double **f = atom->f; int *mask = atom->mask; - imageint *image = atom->image; + tagint *image = atom->image; int nlocal = atom->nlocal; double dx, dy, dz; double unwrap[3]; espring = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->unmap(x[i],image[i],unwrap); dx = unwrap[0] - xoriginal[i][0]; dy = unwrap[1] - xoriginal[i][1]; 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); } espring *= 0.5; } /* ---------------------------------------------------------------------- */ void FixTISpring::post_force_respa(int vflag, int ilevel, int iloop) { if (ilevel == nlevels_respa-1) post_force(vflag); } /* ---------------------------------------------------------------------- */ void FixTISpring::min_post_force(int vflag) { post_force(vflag); } /* ---------------------------------------------------------------------- */ -void FixTISpring::initial_integrate(int vflag) -{ +void FixTISpring::initial_integrate(int vflag) +{ // Update the coupling parameter value. - const bigint t = update->ntimestep - (t0+t_equil); - const double r_switch = 1.0/t_switch; + double t = update->ntimestep - (t0+t_equil); if( (t >= 0) && (t <= t_switch) ) { - lambda = switch_func(t*r_switch); - dlambda = dswitch_func(t*r_switch); + lambda = switch_func(t/t_switch); + dlambda = dswitch_func(t/t_switch); } if( (t >= t_equil+t_switch) && (t <= (t_equil+2*t_switch)) ) { - lambda = switch_func(1.0 - (t - t_switch - t_equil)*r_switch); - dlambda = - dswitch_func(1.0 - (t - t_switch - t_equil)*r_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 ); } -} +} /* ---------------------------------------------------------------------- energy of stretched springs ------------------------------------------------------------------------- */ double FixTISpring::compute_scalar() { double all; MPI_Allreduce(&espring,&all,1,MPI_DOUBLE,MPI_SUM,world); return all; } /* ---------------------------------------------------------------------- information about coupling parameter ------------------------------------------------------------------------- */ double FixTISpring::compute_vector(int n) { linfo[0] = lambda; linfo[1] = dlambda; return linfo[n]; } /* ---------------------------------------------------------------------- memory usage of local atom-based array ------------------------------------------------------------------------- */ double FixTISpring::memory_usage() { double bytes = atom->nmax*3 * sizeof(double); return bytes; } /* ---------------------------------------------------------------------- allocate atom-based array ------------------------------------------------------------------------- */ void FixTISpring::grow_arrays(int nmax) { memory->grow(xoriginal,nmax,3,"fix_ti/spring:xoriginal"); } /* ---------------------------------------------------------------------- copy values within local atom-based array ------------------------------------------------------------------------- */ void FixTISpring::copy_arrays(int i, int j, int delflag) { xoriginal[j][0] = xoriginal[i][0]; xoriginal[j][1] = xoriginal[i][1]; xoriginal[j][2] = xoriginal[i][2]; } /* ---------------------------------------------------------------------- pack values in local atom-based array for exchange with another proc ------------------------------------------------------------------------- */ int FixTISpring::pack_exchange(int i, double *buf) { buf[0] = xoriginal[i][0]; buf[1] = xoriginal[i][1]; buf[2] = xoriginal[i][2]; return 3; } /* ---------------------------------------------------------------------- unpack values in local atom-based array from exchange with another proc ------------------------------------------------------------------------- */ int FixTISpring::unpack_exchange(int nlocal, double *buf) { xoriginal[nlocal][0] = buf[0]; xoriginal[nlocal][1] = buf[1]; xoriginal[nlocal][2] = buf[2]; return 3; } /* ---------------------------------------------------------------------- pack values in local atom-based arrays for restart file ------------------------------------------------------------------------- */ int FixTISpring::pack_restart(int i, double *buf) { buf[0] = 4; buf[1] = xoriginal[i][0]; buf[2] = xoriginal[i][1]; buf[3] = xoriginal[i][2]; return 4; } /* ---------------------------------------------------------------------- unpack values from atom->extra array to restart the fix ------------------------------------------------------------------------- */ void FixTISpring::unpack_restart(int nlocal, int nth) { double **extra = atom->extra; // skip to Nth set of extra values - + int m = 0; for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]); m++; - + xoriginal[nlocal][0] = extra[nlocal][m++]; xoriginal[nlocal][1] = extra[nlocal][m++]; xoriginal[nlocal][2] = extra[nlocal][m++]; } /* ---------------------------------------------------------------------- maxsize of any atom's restart data ------------------------------------------------------------------------- */ int FixTISpring::maxsize_restart() { return 4; } /* ---------------------------------------------------------------------- size of atom nlocal's restart data ------------------------------------------------------------------------- */ int FixTISpring::size_restart(int nlocal) { return 4; } /* ---------------------------------------------------------------------- Switching function. ------------------------------------------------------------------------- */ -double FixTISpring::switch_func(double t) +double FixTISpring::switch_func(double t) { if (sf == 1) return t; 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); } /* ---------------------------------------------------------------------- Switching function derivative. ------------------------------------------------------------------------- */ -double FixTISpring::dswitch_func(double t) +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; } diff --git a/src/USER-MISC/fix_ti_spring.h b/src/USER-MISC/fix_ti_spring.h index bce342162..9b956e707 100755 --- a/src/USER-MISC/fix_ti_spring.h +++ b/src/USER-MISC/fix_ti_spring.h @@ -1,94 +1,95 @@ -/* -*- 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. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing authors: - Rodrigo Freitas (Unicamp/Brazil) - rodrigohb@gmail.com + Contributing authors: + Rodrigo Freitas (UC Berkeley) - rodrigof@berkeley.edu + Mark Asta (UC Berkeley) - mdasta@berkeley.edu Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br ------------------------------------------------------------------------- */ #ifdef FIX_CLASS FixStyle(ti/spring,FixTISpring) #else #ifndef LMP_FIX_TI_SPRING_H #define LMP_FIX_TI_SPRING_H #include "fix.h" namespace LAMMPS_NS { class FixTISpring : public Fix { public: FixTISpring(class LAMMPS *, int, char **); ~FixTISpring(); int setmask(); void init(); void setup(int); void min_setup(int); void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); void initial_integrate(int); double compute_scalar(); double compute_vector(int); - + double memory_usage(); void grow_arrays(int); void copy_arrays(int, int, int); int pack_exchange(int, double *); int unpack_exchange(int, double *); int pack_restart(int, double *); void unpack_restart(int, int); int size_restart(int); int maxsize_restart(); private: double switch_func(double); // Switching function. double dswitch_func(double); // Switching function derivative. double k; // Spring constant. double espring; // Springs energies. double **xoriginal; // Original coords of atoms. double lambda; // Coupling parameter. double dlambda; // Lambda variation with t. double linfo[2]; // Current lambda status. - bigint t_switch; // Total switching steps. - bigint t_equil; // Equilibration time. - bigint t0; // Initial time. + int t_switch; // Total switching steps. + int t_equil; // Equilibration time. + int t0; // Initial time. int sf; // Switching function option. int nlevels_respa; }; } #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: 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. */