forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1607 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
c939c91289
commit
62590294d6
|
|
@ -5,13 +5,9 @@ if ($1 == 1) then
|
||||||
cp style_dipole.h ..
|
cp style_dipole.h ..
|
||||||
|
|
||||||
cp atom_vec_dipole.cpp ..
|
cp atom_vec_dipole.cpp ..
|
||||||
cp compute_temp_dipole.cpp ..
|
|
||||||
cp fix_nve_dipole.cpp ..
|
|
||||||
cp pair_dipole_cut.cpp ..
|
cp pair_dipole_cut.cpp ..
|
||||||
|
|
||||||
cp atom_vec_dipole.h ..
|
cp atom_vec_dipole.h ..
|
||||||
cp compute_temp_dipole.h ..
|
|
||||||
cp fix_nve_dipole.h ..
|
|
||||||
cp pair_dipole_cut.h ..
|
cp pair_dipole_cut.h ..
|
||||||
|
|
||||||
else if ($1 == 0) then
|
else if ($1 == 0) then
|
||||||
|
|
@ -20,13 +16,9 @@ else if ($1 == 0) then
|
||||||
touch ../style_dipole.h
|
touch ../style_dipole.h
|
||||||
|
|
||||||
rm ../atom_vec_dipole.cpp
|
rm ../atom_vec_dipole.cpp
|
||||||
rm ../compute_temp_dipole.cpp
|
|
||||||
rm ../fix_nve_dipole.cpp
|
|
||||||
rm ../pair_dipole_cut.cpp
|
rm ../pair_dipole_cut.cpp
|
||||||
|
|
||||||
rm ../atom_vec_dipole.h
|
rm ../atom_vec_dipole.h
|
||||||
rm ../compute_temp_dipole.h
|
|
||||||
rm ../fix_nve_dipole.h
|
|
||||||
rm ../pair_dipole_cut.h
|
rm ../pair_dipole_cut.h
|
||||||
|
|
||||||
endif
|
endif
|
||||||
|
|
|
||||||
|
|
@ -1,150 +0,0 @@
|
||||||
/* ----------------------------------------------------------------------
|
|
||||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
|
||||||
www.cs.sandia.gov/~sjplimp/lammps.html
|
|
||||||
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
|
|
||||||
|
|
||||||
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.
|
|
||||||
------------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
#include "mpi.h"
|
|
||||||
#include "compute_temp_dipole.h"
|
|
||||||
#include "atom.h"
|
|
||||||
#include "force.h"
|
|
||||||
#include "domain.h"
|
|
||||||
#include "group.h"
|
|
||||||
#include "modify.h"
|
|
||||||
#include "fix.h"
|
|
||||||
#include "error.h"
|
|
||||||
|
|
||||||
using namespace LAMMPS_NS;
|
|
||||||
|
|
||||||
#define INVOKED_SCALAR 1
|
|
||||||
#define INVOKED_VECTOR 2
|
|
||||||
#define INERTIA 0.4 // moment of inertia for a sphere
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
ComputeTempDipole::ComputeTempDipole(LAMMPS *lmp, int narg, char **arg) :
|
|
||||||
Compute(lmp, narg, arg)
|
|
||||||
{
|
|
||||||
if (narg != 3) error->all("Illegal compute temp/dipole command");
|
|
||||||
|
|
||||||
if (!atom->omega_flag || atom->shape == NULL)
|
|
||||||
error->all("Compute temp/dipole requires atom attributes omega, shape");
|
|
||||||
|
|
||||||
scalar_flag = vector_flag = 1;
|
|
||||||
size_vector = 6;
|
|
||||||
extscalar = 0;
|
|
||||||
extvector = 1;
|
|
||||||
tempflag = 1;
|
|
||||||
|
|
||||||
vector = new double[6];
|
|
||||||
inertia = new double[atom->ntypes + 1];
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
ComputeTempDipole::~ComputeTempDipole()
|
|
||||||
{
|
|
||||||
delete [] vector;
|
|
||||||
delete [] inertia;
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
void ComputeTempDipole::init()
|
|
||||||
{
|
|
||||||
fix_dof = 0;
|
|
||||||
for (int i = 0; i < modify->nfix; i++)
|
|
||||||
fix_dof += modify->fix[i]->dof(igroup);
|
|
||||||
recount();
|
|
||||||
|
|
||||||
// moment of inertia for each particle type
|
|
||||||
|
|
||||||
double *mass = atom->mass;
|
|
||||||
double **shape = atom->shape;
|
|
||||||
|
|
||||||
for (int i = 1; i <= atom->ntypes; i++)
|
|
||||||
inertia[i] = INERTIA * mass[i] * 0.25*shape[i][0]*shape[i][0];
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
void ComputeTempDipole::recount()
|
|
||||||
{
|
|
||||||
double natoms = group->count(igroup);
|
|
||||||
dof = 2.0 * domain->dimension * natoms;
|
|
||||||
dof -= extra_dof + fix_dof;
|
|
||||||
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
|
|
||||||
else tfactor = 0.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
double ComputeTempDipole::compute_scalar()
|
|
||||||
{
|
|
||||||
invoked |= INVOKED_SCALAR;
|
|
||||||
|
|
||||||
double **v = atom->v;
|
|
||||||
double *mass = atom->mass;
|
|
||||||
double **omega = atom->omega;
|
|
||||||
int *type = atom->type;
|
|
||||||
int *mask = atom->mask;
|
|
||||||
int nlocal = atom->nlocal;
|
|
||||||
|
|
||||||
// rotational and translational kinetic energy
|
|
||||||
|
|
||||||
double t = 0.0;
|
|
||||||
for (int i = 0; i < nlocal; i++)
|
|
||||||
if (mask[i] & groupbit)
|
|
||||||
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
|
|
||||||
mass[type[i]]
|
|
||||||
+ (omega[i][0] * omega[i][0] + omega[i][1] * omega[i][1] +
|
|
||||||
omega[i][2] * omega[i][2]) * inertia[type[i]];
|
|
||||||
|
|
||||||
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
|
|
||||||
if (dynamic) recount();
|
|
||||||
scalar *= tfactor;
|
|
||||||
return scalar;
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
void ComputeTempDipole::compute_vector()
|
|
||||||
{
|
|
||||||
int i;
|
|
||||||
|
|
||||||
invoked |= INVOKED_VECTOR;
|
|
||||||
|
|
||||||
double **v = atom->v;
|
|
||||||
double *mass = atom->mass;
|
|
||||||
double **omega = atom->omega;
|
|
||||||
int *type = atom->type;
|
|
||||||
int *mask = atom->mask;
|
|
||||||
int nlocal = atom->nlocal;
|
|
||||||
|
|
||||||
double rmass,imass,t[6];
|
|
||||||
for (i = 0; i < 6; i++) t[i] = 0.0;
|
|
||||||
|
|
||||||
// rotational and translational kinetic energy
|
|
||||||
|
|
||||||
for (i = 0; i < nlocal; i++)
|
|
||||||
if (mask[i] & groupbit) {
|
|
||||||
rmass = mass[type[i]];
|
|
||||||
imass = inertia[type[i]];
|
|
||||||
t[0] += rmass*v[i][0]*v[i][0] + imass*omega[i][0]*omega[i][0];
|
|
||||||
t[1] += rmass*v[i][1]*v[i][1] + imass*omega[i][1]*omega[i][1];
|
|
||||||
t[2] += rmass*v[i][2]*v[i][2] + imass*omega[i][2]*omega[i][2];
|
|
||||||
t[3] += rmass*v[i][0]*v[i][1] + imass*omega[i][0]*omega[i][1];
|
|
||||||
t[4] += rmass*v[i][0]*v[i][2] + imass*omega[i][0]*omega[i][2];
|
|
||||||
t[5] += rmass*v[i][1]*v[i][2] + imass*omega[i][1]*omega[i][2];
|
|
||||||
}
|
|
||||||
|
|
||||||
MPI_Allreduce(&t,&vector,6,MPI_DOUBLE,MPI_SUM,world);
|
|
||||||
for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
|
|
||||||
}
|
|
||||||
|
|
@ -1,39 +0,0 @@
|
||||||
/* ----------------------------------------------------------------------
|
|
||||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
|
||||||
www.cs.sandia.gov/~sjplimp/lammps.html
|
|
||||||
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
|
|
||||||
|
|
||||||
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.
|
|
||||||
------------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
#ifndef COMPUTE_TEMP_DIPOLE_H
|
|
||||||
#define COMPUTE_TEMP_DIPOLE_H
|
|
||||||
|
|
||||||
#include "compute.h"
|
|
||||||
|
|
||||||
namespace LAMMPS_NS {
|
|
||||||
|
|
||||||
class ComputeTempDipole : public Compute {
|
|
||||||
public:
|
|
||||||
ComputeTempDipole(class LAMMPS *, int, char **);
|
|
||||||
~ComputeTempDipole();
|
|
||||||
void init();
|
|
||||||
double compute_scalar();
|
|
||||||
void compute_vector();
|
|
||||||
|
|
||||||
private:
|
|
||||||
int fix_dof;
|
|
||||||
double tfactor;
|
|
||||||
double *inertia;
|
|
||||||
|
|
||||||
void recount();
|
|
||||||
};
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
#endif
|
|
||||||
|
|
@ -1,201 +0,0 @@
|
||||||
/* ----------------------------------------------------------------------
|
|
||||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
|
||||||
www.cs.sandia.gov/~sjplimp/lammps.html
|
|
||||||
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
|
|
||||||
|
|
||||||
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.
|
|
||||||
------------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
#include "math.h"
|
|
||||||
#include "string.h"
|
|
||||||
#include "fix_nve_dipole.h"
|
|
||||||
#include "atom.h"
|
|
||||||
#include "atom_vec.h"
|
|
||||||
#include "force.h"
|
|
||||||
#include "update.h"
|
|
||||||
#include "respa.h"
|
|
||||||
#include "error.h"
|
|
||||||
|
|
||||||
using namespace LAMMPS_NS;
|
|
||||||
|
|
||||||
// moment of inertia for a sphere
|
|
||||||
|
|
||||||
#define INERTIA 0.4
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
FixNVEDipole::FixNVEDipole(LAMMPS *lmp, int narg, char **arg) :
|
|
||||||
Fix(lmp, narg, arg)
|
|
||||||
{
|
|
||||||
if (narg != 3) error->all("Illegal fix nve/dipole command");
|
|
||||||
if (!atom->mu_flag || !atom->omega_flag ||
|
|
||||||
!atom->torque_flag || !atom->avec->shape_type)
|
|
||||||
error->all("Fix nve/dipole requires atom attributes "
|
|
||||||
"mu, omega, torque, shape");
|
|
||||||
inertia = new double[atom->ntypes+1];
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
FixNVEDipole::~FixNVEDipole()
|
|
||||||
{
|
|
||||||
delete [] inertia;
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
int FixNVEDipole::setmask()
|
|
||||||
{
|
|
||||||
int mask = 0;
|
|
||||||
mask |= INITIAL_INTEGRATE;
|
|
||||||
mask |= FINAL_INTEGRATE;
|
|
||||||
mask |= INITIAL_INTEGRATE_RESPA;
|
|
||||||
mask |= FINAL_INTEGRATE_RESPA;
|
|
||||||
return mask;
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
void FixNVEDipole::init()
|
|
||||||
{
|
|
||||||
dtv = update->dt;
|
|
||||||
dtf = 0.5 * update->dt * force->ftm2v;
|
|
||||||
|
|
||||||
if (strcmp(update->integrate_style,"respa") == 0)
|
|
||||||
step_respa = ((Respa *) update->integrate)->step;
|
|
||||||
|
|
||||||
// moment of inertia for each particle type
|
|
||||||
|
|
||||||
double *mass = atom->mass;
|
|
||||||
double **shape = atom->shape;
|
|
||||||
|
|
||||||
for (int i = 1; i <= atom->ntypes; i++)
|
|
||||||
inertia[i] = INERTIA * mass[i] * 0.25*shape[i][0]*shape[i][0];
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
void FixNVEDipole::initial_integrate(int vflag)
|
|
||||||
{
|
|
||||||
double dtfm,msq,scale;
|
|
||||||
|
|
||||||
double **x = atom->x;
|
|
||||||
double **v = atom->v;
|
|
||||||
double **f = atom->f;
|
|
||||||
double **mu = atom->mu;
|
|
||||||
double **omega = atom->omega;
|
|
||||||
double **torque = atom->torque;
|
|
||||||
double *mass = atom->mass;
|
|
||||||
double *dipole = atom->dipole;
|
|
||||||
int *type = atom->type;
|
|
||||||
int *mask = atom->mask;
|
|
||||||
int nlocal = atom->nlocal;
|
|
||||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
|
||||||
|
|
||||||
double g[3];
|
|
||||||
|
|
||||||
// update v,x for all particles
|
|
||||||
// update omega,mu for all dipoles
|
|
||||||
// d_omega/dt = torque / inertia
|
|
||||||
// d_mu/dt = omega cross mu
|
|
||||||
// renormalize mu to dipole length
|
|
||||||
|
|
||||||
for (int i = 0; i < nlocal; i++) {
|
|
||||||
if (mask[i] & groupbit) {
|
|
||||||
dtfm = dtf / mass[type[i]];
|
|
||||||
v[i][0] += dtfm * f[i][0];
|
|
||||||
v[i][1] += dtfm * f[i][1];
|
|
||||||
v[i][2] += dtfm * f[i][2];
|
|
||||||
x[i][0] += dtv * v[i][0];
|
|
||||||
x[i][1] += dtv * v[i][1];
|
|
||||||
x[i][2] += dtv * v[i][2];
|
|
||||||
if (dipole[type[i]] > 0.0) {
|
|
||||||
dtfm = dtf / inertia[type[i]];
|
|
||||||
omega[i][0] += dtfm * torque[i][0];
|
|
||||||
omega[i][1] += dtfm * torque[i][1];
|
|
||||||
omega[i][2] += dtfm * torque[i][2];
|
|
||||||
|
|
||||||
g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2] - omega[i][2]*mu[i][1]);
|
|
||||||
g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0] - omega[i][0]*mu[i][2]);
|
|
||||||
g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1] - omega[i][1]*mu[i][0]);
|
|
||||||
|
|
||||||
msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
|
|
||||||
scale = dipole[type[i]]/sqrt(msq);
|
|
||||||
mu[i][0] = g[0]*scale;
|
|
||||||
mu[i][1] = g[1]*scale;
|
|
||||||
mu[i][2] = g[2]*scale;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
void FixNVEDipole::final_integrate()
|
|
||||||
{
|
|
||||||
double dtfm;
|
|
||||||
|
|
||||||
double **v = atom->v;
|
|
||||||
double **f = atom->f;
|
|
||||||
double **omega = atom->omega;
|
|
||||||
double **torque = atom->torque;
|
|
||||||
double *mass = atom->mass;
|
|
||||||
double *dipole = atom->dipole;
|
|
||||||
int *type = atom->type;
|
|
||||||
int *mask = atom->mask;
|
|
||||||
int nlocal = atom->nlocal;
|
|
||||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
|
||||||
|
|
||||||
// update v for all particles
|
|
||||||
// update omega for all dipoles
|
|
||||||
|
|
||||||
for (int i = 0; i < nlocal; i++) {
|
|
||||||
if (mask[i] & groupbit) {
|
|
||||||
dtfm = dtf / mass[type[i]];
|
|
||||||
v[i][0] += dtfm * f[i][0];
|
|
||||||
v[i][1] += dtfm * f[i][1];
|
|
||||||
v[i][2] += dtfm * f[i][2];
|
|
||||||
if (dipole[type[i]] > 0.0) {
|
|
||||||
dtfm = dtf / inertia[type[i]];
|
|
||||||
omega[i][0] += dtfm * torque[i][0];
|
|
||||||
omega[i][1] += dtfm * torque[i][1];
|
|
||||||
omega[i][2] += dtfm * torque[i][2];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
void FixNVEDipole::initial_integrate_respa(int vflag, int ilevel, int flag)
|
|
||||||
{
|
|
||||||
if (flag) return; // only used by NPT,NPH
|
|
||||||
|
|
||||||
dtv = step_respa[ilevel];
|
|
||||||
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
|
|
||||||
|
|
||||||
if (ilevel == 0) initial_integrate(vflag);
|
|
||||||
else final_integrate();
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
void FixNVEDipole::final_integrate_respa(int ilevel)
|
|
||||||
{
|
|
||||||
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
|
|
||||||
final_integrate();
|
|
||||||
}
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
void FixNVEDipole::reset_dt()
|
|
||||||
{
|
|
||||||
dtv = update->dt;
|
|
||||||
dtf = 0.5 * update->dt * force->ftm2v;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
@ -1,41 +0,0 @@
|
||||||
/* ----------------------------------------------------------------------
|
|
||||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
|
||||||
www.cs.sandia.gov/~sjplimp/lammps.html
|
|
||||||
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
|
|
||||||
|
|
||||||
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.
|
|
||||||
------------------------------------------------------------------------- */
|
|
||||||
|
|
||||||
#ifndef FIX_NVE_DIPOLE_H
|
|
||||||
#define FIX_NVE_DIPOLE_H
|
|
||||||
|
|
||||||
#include "fix.h"
|
|
||||||
|
|
||||||
namespace LAMMPS_NS {
|
|
||||||
|
|
||||||
class FixNVEDipole : public Fix {
|
|
||||||
public:
|
|
||||||
FixNVEDipole(class LAMMPS *, int, char **);
|
|
||||||
~FixNVEDipole();
|
|
||||||
int setmask();
|
|
||||||
void init();
|
|
||||||
void initial_integrate(int);
|
|
||||||
void final_integrate();
|
|
||||||
void initial_integrate_respa(int, int, int);
|
|
||||||
void final_integrate_respa(int);
|
|
||||||
void reset_dt();
|
|
||||||
|
|
||||||
private:
|
|
||||||
double dtv,dtf;
|
|
||||||
double *step_respa;
|
|
||||||
double *inertia;
|
|
||||||
};
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
#endif
|
|
||||||
|
|
@ -19,22 +19,6 @@
|
||||||
AtomStyle(dipole,AtomVecDipole)
|
AtomStyle(dipole,AtomVecDipole)
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
#ifdef ComputeInclude
|
|
||||||
#include "compute_temp_dipole.h"
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#ifdef ComputeClass
|
|
||||||
ComputeStyle(temp/dipole,ComputeTempDipole)
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#ifdef FixInclude
|
|
||||||
#include "fix_nve_dipole.h"
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#ifdef FixClass
|
|
||||||
FixStyle(nve/dipole,FixNVEDipole)
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#ifdef PairInclude
|
#ifdef PairInclude
|
||||||
#include "pair_dipole_cut.h"
|
#include "pair_dipole_cut.h"
|
||||||
#endif
|
#endif
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue