From 05649f118942a05cfede3359b74a43ed28d99ad5 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 7 Nov 2013 14:38:21 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10987 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-MISC/fix_ti_rs.cpp | 394 ++++++++--------- src/USER-MISC/fix_ti_rs.h | 132 +++--- src/USER-MISC/fix_ti_spring.cpp | 730 ++++++++++++++++---------------- src/USER-MISC/fix_ti_spring.h | 188 ++++---- 4 files changed, 722 insertions(+), 722 deletions(-) diff --git a/src/USER-MISC/fix_ti_rs.cpp b/src/USER-MISC/fix_ti_rs.cpp index e4e78fcd1f..f80fc4cb17 100755 --- a/src/USER-MISC/fix_ti_rs.cpp +++ b/src/USER-MISC/fix_ti_rs.cpp @@ -1,197 +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 = atoi(arg[5]); - t_equil = atoi(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; -} +/* ------------------------------------------------------------------------- + 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 = atoi(arg[5]); + t_equil = atoi(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 index fabe1b49ee..73aac11aed 100755 --- a/src/USER-MISC/fix_ti_rs.h +++ b/src/USER-MISC/fix_ti_rs.h @@ -1,66 +1,66 @@ -/* ---------------------------------------------------------------------- - 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. - 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 +/* ---------------------------------------------------------------------- + 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. + 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 index 3fe1715e50..a2ddb9a0f0 100755 --- a/src/USER-MISC/fix_ti_spring.cpp +++ b/src/USER-MISC/fix_ti_spring.cpp @@ -1,365 +1,365 @@ -/* ------------------------------------------------------------------------- - 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" -#include "force.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; - 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; - 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 = atoi(arg[4]); // Switching time. - t_equil = atoi(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); - - 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. - int t = update->ntimestep - t0; - 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]; - - 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) -{ - // 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, 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 (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; -} +/* ------------------------------------------------------------------------- + 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" +#include "force.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; + 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; + 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 = atoi(arg[4]); // Switching time. + t_equil = atoi(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); + + 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. + int t = update->ntimestep - t0; + 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]; + + 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) +{ + // 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, 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 (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 index 41b2b91b0f..76c7479bb0 100755 --- a/src/USER-MISC/fix_ti_spring.h +++ b/src/USER-MISC/fix_ti_spring.h @@ -1,94 +1,94 @@ -/* ---------------------------------------------------------------------- - 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/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. - 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. - -*/ +/* ---------------------------------------------------------------------- + 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/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. + 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. + +*/