diff --git a/src/USER-MISC/README b/src/USER-MISC/README index ba19f77ac0546ff172955cf4261cf48e992049ec..8e83a4414f14e826d1f0d5a87d7c84845d9107af 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -37,6 +37,8 @@ dihedral_style table, Andrew Jewett, jewett.aij@gmail.com, 10 Jan 12 fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 fix imd, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 2009 fix smd, Axel Kohlmeyer, akohlmey at gmail.com, 19 May 2008 +fix ti/rs, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013 +fix ti/spring, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013 improper_style cossq, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12 improper_style fourier, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 improper_style ring, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12 diff --git a/src/USER-MISC/fix_ti_rs.cpp b/src/USER-MISC/fix_ti_rs.cpp new file mode 100755 index 0000000000000000000000000000000000000000..8ca845c87f24c48a6a31b2d17b02b532f8a4e779 --- /dev/null +++ b/src/USER-MISC/fix_ti_rs.cpp @@ -0,0 +1,197 @@ +/* ------------------------------------------------------------------------- + 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" + +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 = atof(arg[5]); + t_equil = atof(arg[6]); + t0 = update->ntimestep; + 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"); + + // Coupling parameter limits and initialization. + l_initial = atof(arg[3]); + l_final = atof(arg[4]); + sf = 1; + if (narg > 7) { + if (strcmp(arg[7], "function") == 0) sf = atoi(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. + 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); + } + + 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 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 new file mode 100755 index 0000000000000000000000000000000000000000..f0c8bdbe904d86f4989d77d2ebcefff6b6b2418b --- /dev/null +++ b/src/USER-MISC/fix_ti_rs.h @@ -0,0 +1,60 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(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. + 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 diff --git a/src/USER-MISC/fix_ti_spring.cpp b/src/USER-MISC/fix_ti_spring.cpp new file mode 100755 index 0000000000000000000000000000000000000000..ecf76fba8a6b6084d0aac0291f3a6aa2ee7a6a29 --- /dev/null +++ b/src/USER-MISC/fix_ti_spring.cpp @@ -0,0 +1,359 @@ +/* ------------------------------------------------------------------------- + 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 "fix_ti_spring.h" +#include "atom.h" +#include "update.h" +#include "domain.h" +#include "respa.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +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; + vector_flag = 1; + size_vector = 2; + global_freq = 1; + extscalar = 1; + extvector = 1; + + // Spring constant. + k = atof(arg[3]); + espring = 0.0; + if (k <= 0.0) error->all(FLERR,"Illegal fix ti/spring command"); + + // Initial position. + xoriginal = NULL; + grow_arrays(atom->nmax); + atom->add_callback(0); + atom->add_callback(1); + + double **x = atom->x; + int *mask = atom->mask; + 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 = atof(arg[4]); // Switching time. + t_equil = atof(arg[5]); // Equilibration time. + 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"); + + // Coupling parameter initialization. + sf = 1; + if (narg > 6) { + 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)) + error->all(FLERR,"Illegal fix ti/spring switching function"); + } + lambda = switch_func(0); + dlambda = dswitch_func(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. + int t = update->ntimestep - t0; + espring = 0.0; + if(t < t_equil) return; + + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + tagint *image = atom->image; + int nlocal = atom->nlocal; + + double dx, dy, dz; + double unwrap[3]; + + 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) +{ + // Update the coupling parameter value. + 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); + } + + 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 ); + } +} + +/* ---------------------------------------------------------------------- + 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) +{ + 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) +{ + 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) +{ + 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 new file mode 100755 index 0000000000000000000000000000000000000000..21d7ef5d14468126c3dc1a00947840c0e891e281 --- /dev/null +++ b/src/USER-MISC/fix_ti_spring.h @@ -0,0 +1,88 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(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 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. + 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. + +*/