Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81167077
fix_ti_spring.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Sep 4, 13:14
Size
10 KB
Mime Type
text/x-c
Expires
Fri, Sep 6, 13:14 (2 d)
Engine
blob
Format
Raw Data
Handle
18793668
Attached To
rLAMMPS lammps
fix_ti_spring.cpp
View Options
/* -------------------------------------------------------------------------
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;
imageint *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.
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");
// Coupling parameter initialization.
sf = 1;
if (narg > 6) {
if (strcmp(arg[6], "function") == 0) sf = force->inumeric(FLERR,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.
bigint 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;
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.
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);
}
}
/* ----------------------------------------------------------------------
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)
{
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;
}
Event Timeline
Log In to Comment