From 116481677601d0bdc148cad1be5f4b58da7c004d Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 20 Jun 2007 13:32:19 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@624 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/ASPHERE/fix_npt_asphere.cpp | 300 ++++++++ .../fix_npt_asphere.h} | 29 +- src/ASPHERE/fix_nvt_asphere.cpp | 237 ++++++ src/ASPHERE/fix_nvt_asphere.h | 39 + src/DIPOLE/Install.csh | 32 + src/DIPOLE/atom_vec_dipole.cpp | 657 +++++++++++++++++ src/DIPOLE/atom_vec_dipole.h | 58 ++ src/DIPOLE/compute_temp_dipole.cpp | 144 ++++ src/DIPOLE/compute_temp_dipole.h | 39 + src/DIPOLE/fix_nve_dipole.cpp | 189 +++++ src/DIPOLE/fix_nve_dipole.h | 40 ++ src/DIPOLE/pair_dipole_cut.cpp | 625 ++++++++++++++++ src/DIPOLE/pair_dipole_cut.h | 48 ++ src/DIPOLE/style_dipole.h | 44 ++ src/MANYBODY/Install.csh | 4 - src/MANYBODY/style_manybody.h | 2 - src/compute_temp_deform.cpp | 193 +++++ src/compute_temp_deform.h | 38 + src/fix_deform.cpp | 674 ++++++++++++++++++ src/{fix_uniaxial.h => fix_deform.h} | 53 +- src/fix_nvt_sllod.cpp | 278 ++++++++ src/fix_nvt_sllod.h | 33 + src/fix_uniaxial.cpp | 212 ------ src/fix_volume_rescale.cpp | 194 ----- src/math_extra.h | 435 +++++++++++ src/style_asphere.h | 0 src/style_class2.h | 56 -- src/style_colloid.h | 0 src/style_dipole.h | 0 src/style_dpd.h | 28 - src/style_granular.h | 50 -- src/style_kspace.h | 38 - src/style_manybody.h | 30 - src/style_molecule.h | 116 --- src/style_opt.h | 30 - src/style_xtc.h | 20 - 36 files changed, 4150 insertions(+), 815 deletions(-) create mode 100755 src/ASPHERE/fix_npt_asphere.cpp rename src/{fix_volume_rescale.h => ASPHERE/fix_npt_asphere.h} (55%) mode change 100644 => 100755 create mode 100755 src/ASPHERE/fix_nvt_asphere.cpp create mode 100755 src/ASPHERE/fix_nvt_asphere.h create mode 100644 src/DIPOLE/Install.csh create mode 100644 src/DIPOLE/atom_vec_dipole.cpp create mode 100644 src/DIPOLE/atom_vec_dipole.h create mode 100644 src/DIPOLE/compute_temp_dipole.cpp create mode 100644 src/DIPOLE/compute_temp_dipole.h create mode 100644 src/DIPOLE/fix_nve_dipole.cpp create mode 100644 src/DIPOLE/fix_nve_dipole.h create mode 100644 src/DIPOLE/pair_dipole_cut.cpp create mode 100644 src/DIPOLE/pair_dipole_cut.h create mode 100644 src/DIPOLE/style_dipole.h create mode 100644 src/compute_temp_deform.cpp create mode 100644 src/compute_temp_deform.h create mode 100644 src/fix_deform.cpp rename src/{fix_uniaxial.h => fix_deform.h} (53%) create mode 100644 src/fix_nvt_sllod.cpp create mode 100644 src/fix_nvt_sllod.h delete mode 100644 src/fix_uniaxial.cpp delete mode 100644 src/fix_volume_rescale.cpp create mode 100755 src/math_extra.h create mode 100644 src/style_asphere.h create mode 100644 src/style_colloid.h create mode 100644 src/style_dipole.h diff --git a/src/ASPHERE/fix_npt_asphere.cpp b/src/ASPHERE/fix_npt_asphere.cpp new file mode 100755 index 0000000000..147fec0cca --- /dev/null +++ b/src/ASPHERE/fix_npt_asphere.cpp @@ -0,0 +1,300 @@ +/* ---------------------------------------------------------------------- + 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 author: Mike Brown (SNL) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "stdlib.h" +#include "math.h" +#include "fix_npt_asphere.h" +#include "math_extra.h" +#include "atom.h" +#include "force.h" +#include "compute.h" +#include "kspace.h" +#include "update.h" +#include "domain.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) : + FixNPT(lmp, narg, arg) +{ + if (atom->quat == NULL || atom->angmom == NULL || atom->torque == NULL) + error->all("Fix npt/asphere requires atom attributes " + "quat, angmom, torque"); +} + +/* ---------------------------------------------------------------------- */ + +void FixNPTASphere::init() +{ + FixNPT::init(); + dtq = 0.5 * update->dt; +} + +/* ---------------------------------------------------------------------- + 1st half of Verlet update +------------------------------------------------------------------------- */ + +void FixNPTASphere::initial_integrate() +{ + int i; + + double delta = update->ntimestep - update->beginstep; + delta /= update->endstep - update->beginstep; + + // update eta_dot + + t_target = t_start + delta * (t_stop-t_start); + f_eta = t_freq*t_freq * (t_current/t_target - 1.0); + eta_dot += f_eta*dthalf; + eta_dot *= drag_factor; + eta += dtv*eta_dot; + + // update omega_dot + // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change + + double f_omega; + double denskt = (atom->natoms*boltz*t_target) / + (domain->xprd*domain->yprd*domain->zprd) * nktv2p; + + for (i = 0; i < 3; i++) { + p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); + f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt; + omega_dot[i] += f_omega*dthalf; + omega_dot[i] *= drag_factor; + omega[i] += dtv*omega_dot[i]; + factor[i] = exp(-dthalf*(eta_dot+omega_dot[i])); + dilation[i] = exp(dthalf*omega_dot[i]); + } + ang_factor = exp(-dthalf*eta_dot); + + // v update only for atoms in NPT group + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double **quat = atom->quat; + double **angmom = atom->angmom; + double **torque = atom->torque; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double dtfm; + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } + } + + // rescale simulation box and all owned atoms by 1/2 step + + box_dilate(0); + + // x update by full step only for atoms in NPT group + + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } + + // update angular momentum by 1/2 step + // update quaternion a full step via Richardson iteration + // returns new normalized quaternion + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + angmom[i][0] = angmom[i][0] * ang_factor + dtf * torque[i][0]; + angmom[i][1] = angmom[i][1] * ang_factor + dtf * torque[i][1]; + angmom[i][2] = angmom[i][2] * ang_factor + dtf * torque[i][2]; + + double inertia[3]; + calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia); + richardson(quat[i],angmom[i],inertia); + } + } + + // rescale simulation box and all owned atoms by 1/2 step + // redo KSpace coeffs since volume has changed + + box_dilate(0); + if (kspace_flag) force->kspace->setup(); +} + +/* ---------------------------------------------------------------------- + 2nd half of Verlet update +------------------------------------------------------------------------- */ + +void FixNPTASphere::final_integrate() +{ + int i; + + // v update only for atoms in NPT group + + double **v = atom->v; + double **f = atom->f; + double **angmom = atom->angmom; + double **torque = atom->torque; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double dtfm; + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + + angmom[i][0] = (angmom[i][0] + dtf * torque[i][0]) * ang_factor; + angmom[i][1] = (angmom[i][1] + dtf * torque[i][1]) * ang_factor; + angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * ang_factor; + } + } + + // compute new T,P + + t_current = temperature->compute_scalar(); + if (press_couple == 0) { + if (ptemperature) double ptmp = ptemperature->compute_scalar(); + double tmp = pressure->compute_scalar(); + } else { + temperature->compute_vector(); + if (ptemperature) ptemperature->compute_vector(); + pressure->compute_vector(); + } + couple(); + + // update eta_dot + + f_eta = t_freq*t_freq * (t_current/t_target - 1.0); + eta_dot += f_eta*dthalf; + eta_dot *= drag_factor; + + // update omega_dot + // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change + + double f_omega; + double denskt = (atom->natoms*boltz*t_target) / + (domain->xprd*domain->yprd*domain->zprd) * nktv2p; + + for (i = 0; i < 3; i++) { + f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt; + omega_dot[i] += f_omega*dthalf; + omega_dot[i] *= drag_factor; + } +} + +/* ---------------------------------------------------------------------- + Richardson iteration to update quaternion accurately +------------------------------------------------------------------------- */ + +void FixNPTASphere::richardson(double *q, double *m, double *moments) +{ + // compute omega at 1/2 step from m at 1/2 step and q at 0 + + double w[3]; + omega_from_mq(q,m,moments,w); + + // full update from dq/dt = 1/2 w q + + double wq[4]; + MathExtra::multiply_vec_quat(w,q,wq); + + double qfull[4]; + qfull[0] = q[0] + dtq * wq[0]; + qfull[1] = q[1] + dtq * wq[1]; + qfull[2] = q[2] + dtq * wq[2]; + qfull[3] = q[3] + dtq * wq[3]; + MathExtra::normalize4(qfull); + + // 1st half of update from dq/dt = 1/2 w q + + double qhalf[4]; + qhalf[0] = q[0] + 0.5*dtq * wq[0]; + qhalf[1] = q[1] + 0.5*dtq * wq[1]; + qhalf[2] = q[2] + 0.5*dtq * wq[2]; + qhalf[3] = q[3] + 0.5*dtq * wq[3]; + MathExtra::normalize4(qhalf); + + // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step + // recompute wq + + omega_from_mq(qhalf,m,moments,w); + MathExtra::multiply_vec_quat(w,qhalf,wq); + + // 2nd half of update from dq/dt = 1/2 w q + + qhalf[0] += 0.5*dtq * wq[0]; + qhalf[1] += 0.5*dtq * wq[1]; + qhalf[2] += 0.5*dtq * wq[2]; + qhalf[3] += 0.5*dtq * wq[3]; + MathExtra::normalize4(qhalf); + + // corrected Richardson update + + q[0] = 2.0*qhalf[0] - qfull[0]; + q[1] = 2.0*qhalf[1] - qfull[1]; + q[2] = 2.0*qhalf[2] - qfull[2]; + q[3] = 2.0*qhalf[3] - qfull[3]; + MathExtra::normalize4(q); +} + +/* ---------------------------------------------------------------------- + compute omega from angular momentum + w = omega = angular velocity in space frame + wbody = angular velocity in body frame + project space-frame angular momentum onto body axes + and divide by principal moments +------------------------------------------------------------------------- */ + +void FixNPTASphere::omega_from_mq(double *q, double *m, double *inertia, + double *w) +{ + double rot[3][3]; + MathExtra::quat_to_mat(q,rot); + + double wbody[3]; + MathExtra::transpose_times_column3(rot,m,wbody); + wbody[0] /= inertia[0]; + wbody[1] /= inertia[1]; + wbody[2] /= inertia[2]; + MathExtra::times_column3(rot,wbody,w); +} + +/* ---------------------------------------------------------------------- + calculate the moment of inertia for an ellipsoid, from mass and radii + shape = x,y,z radii in body frame +------------------------------------------------------------------------- */ + +void FixNPTASphere::calculate_inertia(double mass, double *shape, + double *inertia) +{ + inertia[0] = mass*(shape[1]*shape[1]+shape[2]*shape[2])/5.0; + inertia[1] = mass*(shape[0]*shape[0]+shape[2]*shape[2])/5.0; + inertia[2] = mass*(shape[0]*shape[0]+shape[1]*shape[1])/5.0; +} diff --git a/src/fix_volume_rescale.h b/src/ASPHERE/fix_npt_asphere.h old mode 100644 new mode 100755 similarity index 55% rename from src/fix_volume_rescale.h rename to src/ASPHERE/fix_npt_asphere.h index d19fe08074..83a4da2eb1 --- a/src/fix_volume_rescale.h +++ b/src/ASPHERE/fix_npt_asphere.h @@ -11,29 +11,28 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef FIX_VOL_RESCALE_H -#define FIX_VOL_RESCALE_H +#ifndef FIX_NPT_ASPHERE_H +#define FIX_NPT_ASPHERE_H -#include "fix.h" +#include "fix_npt.h" namespace LAMMPS_NS { -class FixVolRescale : public Fix { +class FixNPTASphere : public FixNPT { public: - FixVolRescale(class LAMMPS *, int, char **); - ~FixVolRescale(); - int setmask(); + FixNPTASphere(class LAMMPS *, int, char **); + ~FixNPTASphere() {} void init(); - void end_of_step(); + void initial_integrate(); + void final_integrate(); private: - int xflag,yflag,zflag; - double xlo_start,xlo_stop,xhi_start,xhi_stop; - double ylo_start,ylo_stop,yhi_start,yhi_stop; - double zlo_start,zlo_stop,zhi_start,zhi_stop; - int kspace_flag; // 1 if KSpace invoked, 0 if not - int nrigid; // number of rigid fixes - int *rfix; // indices of rigid fixes + double dtq; + double ang_factor; + + void richardson(double *, double *, double *); + void omega_from_mq(double *, double *, double *, double *); + void calculate_inertia(double mass, double *shape, double *inertia); }; } diff --git a/src/ASPHERE/fix_nvt_asphere.cpp b/src/ASPHERE/fix_nvt_asphere.cpp new file mode 100755 index 0000000000..df2e8b3bc3 --- /dev/null +++ b/src/ASPHERE/fix_nvt_asphere.cpp @@ -0,0 +1,237 @@ +/* ---------------------------------------------------------------------- + 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 author: Mike Brown (SNL) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "stdlib.h" +#include "math.h" +#include "fix_nvt_asphere.h" +#include "math_extra.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "group.h" +#include "update.h" +#include "modify.h" +#include "compute.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixNVTASphere::FixNVTASphere(LAMMPS *lmp, int narg, char **arg) : + FixNVT(lmp, narg, arg) +{ + if (atom->quat == NULL || atom->angmom == NULL || atom->torque == NULL) + error->all("Fix nvt/asphere requires atom attributes " + "quat, angmom, torque"); +} + +/* ---------------------------------------------------------------------- */ + +void FixNVTASphere::init() +{ + FixNVT::init(); + dtq = 0.5 * update->dt; +} + +/* ---------------------------------------------------------------------- */ + +void FixNVTASphere::initial_integrate() +{ + double dtfm; + + double delta = update->ntimestep - update->beginstep; + delta /= update->endstep - update->beginstep; + t_target = t_start + delta * (t_stop-t_start); + + // update eta_dot + + f_eta = t_freq*t_freq * (t_current/t_target - 1.0); + eta_dot += f_eta*dthalf; + eta_dot *= drag_factor; + eta += dtv*eta_dot; + factor = exp(-dthalf*eta_dot); + + // update v and x of only atoms in NVT group + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double **quat = atom->quat; + double **angmom = atom->angmom; + double **torque = atom->torque; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + 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]; + + // update angular momentum by 1/2 step + // update quaternion a full step via Richardson iteration + // returns new normalized quaternion + + angmom[i][0] = angmom[i][0] * factor + dtf * torque[i][0]; + angmom[i][1] = angmom[i][1] * factor + dtf * torque[i][1]; + angmom[i][2] = angmom[i][2] * factor + dtf * torque[i][2]; + + double inertia[3]; + calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia); + richardson(quat[i],angmom[i],inertia); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNVTASphere::final_integrate() +{ + double dtfm; + + // update v of only atoms in NVT group + + double **v = atom->v; + double **f = atom->f; + double **angmom = atom->angmom; + double **torque = atom->torque; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + + angmom[i][0] = angmom[i][0] * factor + dtf * torque[i][0]; + angmom[i][1] = angmom[i][1] * factor + dtf * torque[i][1]; + angmom[i][2] = angmom[i][2] * factor + dtf * torque[i][2]; + } + } + + // compute current T + + t_current = temperature->compute_scalar(); + + // update eta_dot + + f_eta = t_freq*t_freq * (t_current/t_target - 1.0); + eta_dot += f_eta*dthalf; + eta_dot *= drag_factor; +} + +/* ---------------------------------------------------------------------- + Richardson iteration to update quaternion accurately +------------------------------------------------------------------------- */ + +void FixNVTASphere::richardson(double *q, double *m, double *moments) +{ + // compute omega at 1/2 step from m at 1/2 step and q at 0 + + double w[3]; + omega_from_mq(q,m,moments,w); + + // full update from dq/dt = 1/2 w q + + double wq[4]; + MathExtra::multiply_vec_quat(w,q,wq); + + double qfull[4]; + qfull[0] = q[0] + dtq * wq[0]; + qfull[1] = q[1] + dtq * wq[1]; + qfull[2] = q[2] + dtq * wq[2]; + qfull[3] = q[3] + dtq * wq[3]; + MathExtra::normalize4(qfull); + + // 1st half of update from dq/dt = 1/2 w q + + double qhalf[4]; + qhalf[0] = q[0] + 0.5*dtq * wq[0]; + qhalf[1] = q[1] + 0.5*dtq * wq[1]; + qhalf[2] = q[2] + 0.5*dtq * wq[2]; + qhalf[3] = q[3] + 0.5*dtq * wq[3]; + MathExtra::normalize4(qhalf); + + // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step + // recompute wq + + omega_from_mq(qhalf,m,moments,w); + MathExtra::multiply_vec_quat(w,qhalf,wq); + + // 2nd half of update from dq/dt = 1/2 w q + + qhalf[0] += 0.5*dtq * wq[0]; + qhalf[1] += 0.5*dtq * wq[1]; + qhalf[2] += 0.5*dtq * wq[2]; + qhalf[3] += 0.5*dtq * wq[3]; + MathExtra::normalize4(qhalf); + + // corrected Richardson update + + q[0] = 2.0*qhalf[0] - qfull[0]; + q[1] = 2.0*qhalf[1] - qfull[1]; + q[2] = 2.0*qhalf[2] - qfull[2]; + q[3] = 2.0*qhalf[3] - qfull[3]; + MathExtra::normalize4(q); +} + +/* ---------------------------------------------------------------------- + compute omega from angular momentum + w = omega = angular velocity in space frame + wbody = angular velocity in body frame + project space-frame angular momentum onto body axes + and divide by principal moments +------------------------------------------------------------------------- */ + +void FixNVTASphere::omega_from_mq(double *q, double *m, double *inertia, + double *w) +{ + double rot[3][3]; + MathExtra::quat_to_mat(q,rot); + + double wbody[3]; + MathExtra::transpose_times_column3(rot,m,wbody); + wbody[0] /= inertia[0]; + wbody[1] /= inertia[1]; + wbody[2] /= inertia[2]; + MathExtra::times_column3(rot,wbody,w); +} + +/* ---------------------------------------------------------------------- + calculate the moment of inertia for an ELLIPSOID, from mass and radii + shape = x,y,z radii in body frame +------------------------------------------------------------------------- */ + +void FixNVTASphere::calculate_inertia(double mass, double *shape, + double *inertia) +{ + inertia[0] = mass*(shape[1]*shape[1]+shape[2]*shape[2])/5.0; + inertia[1] = mass*(shape[0]*shape[0]+shape[2]*shape[2])/5.0; + inertia[2] = mass*(shape[0]*shape[0]+shape[1]*shape[1])/5.0; +} diff --git a/src/ASPHERE/fix_nvt_asphere.h b/src/ASPHERE/fix_nvt_asphere.h new file mode 100755 index 0000000000..d07e5b2d34 --- /dev/null +++ b/src/ASPHERE/fix_nvt_asphere.h @@ -0,0 +1,39 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef FIX_NVT_ASPHERE_H +#define FIX_NVT_ASPHERE_H + +#include "fix_nvt.h" + +namespace LAMMPS_NS { + +class FixNVTASphere : public FixNVT { + public: + FixNVTASphere(class LAMMPS *, int, char **); + ~FixNVTASphere() {} + void init(); + void initial_integrate(); + void final_integrate(); + + private: + double dtq; + + void richardson(double *, double *, double *); + void omega_from_mq(double *, double *, double *, double *); + void calculate_inertia(double mass, double *shape, double *inertia); +}; + +} + +#endif diff --git a/src/DIPOLE/Install.csh b/src/DIPOLE/Install.csh new file mode 100644 index 0000000000..deaadef27c --- /dev/null +++ b/src/DIPOLE/Install.csh @@ -0,0 +1,32 @@ +# Install/unInstall package classes in LAMMPS + +if ($1 == 1) then + + cp style_dipole.h .. + + cp atom_vec_dipole.cpp .. + cp compute_temp_dipole.cpp .. + cp fix_nve_dipole.cpp .. + cp pair_dipole_cut.cpp .. + + cp atom_vec_dipole.h .. + cp compute_temp_dipole.h .. + cp fix_nve_dipole.h .. + cp pair_dipole_cut.h .. + +else if ($1 == 0) then + + rm ../style_dipole.h + touch ../style_dipole.h + + rm ../atom_vec_dipole.cpp + rm ../compute_temp_dipole.cpp + rm ../fix_nve_dipole.cpp + rm ../pair_dipole_cut.cpp + + rm ../atom_vec_dipole.h + rm ../compute_temp_dipole.h + rm ../fix_nve_dipole.h + rm ../pair_dipole_cut.h + +endif diff --git a/src/DIPOLE/atom_vec_dipole.cpp b/src/DIPOLE/atom_vec_dipole.cpp new file mode 100644 index 0000000000..2114bff212 --- /dev/null +++ b/src/DIPOLE/atom_vec_dipole.cpp @@ -0,0 +1,657 @@ +/* ---------------------------------------------------------------------- + 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 "stdlib.h" +#include "atom_vec_dipole.h" +#include "atom.h" +#include "domain.h" +#include "modify.h" +#include "fix.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define DELTA 10000 + +/* ---------------------------------------------------------------------- */ + +AtomVecDipole::AtomVecDipole(LAMMPS *lmp, int narg, char **arg) : + AtomVec(lmp, narg, arg) +{ + mass_type = 1; + shape_type = 1; + dipole_type = 1; + comm_x_only = comm_f_only = 0; + size_comm = 6; + size_reverse = 6; + size_border = 10; + size_data_atom = 9; + size_data_vel = 7; + xcol_data = 4; +} + +/* ---------------------------------------------------------------------- + grow atom arrays + n = 0 grows arrays by DELTA + n > 0 allocates arrays to size n +------------------------------------------------------------------------- */ + +void AtomVecDipole::grow(int n) +{ + if (n == 0) nmax += DELTA; + else nmax = n; + atom->nmax = nmax; + + tag = atom->tag = (int *) + memory->srealloc(atom->tag,nmax*sizeof(int),"atom:tag"); + type = atom->type = (int *) + memory->srealloc(atom->type,nmax*sizeof(int),"atom:type"); + mask = atom->mask = (int *) + memory->srealloc(atom->mask,nmax*sizeof(int),"atom:mask"); + image = atom->image = (int *) + memory->srealloc(atom->image,nmax*sizeof(int),"atom:image"); + x = atom->x = memory->grow_2d_double_array(atom->x,nmax,3,"atom:x"); + v = atom->v = memory->grow_2d_double_array(atom->v,nmax,3,"atom:v"); + f = atom->f = memory->grow_2d_double_array(atom->f,nmax,3,"atom:f"); + + q = atom->q = (double *) + memory->srealloc(atom->q,nmax*sizeof(double),"atom:q"); + mu = atom->mu = + memory->grow_2d_double_array(atom->mu,nmax,3,"atom:mu"); + omega = atom->omega = + memory->grow_2d_double_array(atom->omega,nmax,3,"atom:omega"); + torque = atom->torque = + memory->grow_2d_double_array(atom->torque,nmax,3,"atom:torque"); + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecDipole::copy(int i, int j) +{ + tag[j] = tag[i]; + type[j] = type[i]; + mask[j] = mask[i]; + image[j] = image[i]; + x[j][0] = x[i][0]; + x[j][1] = x[i][1]; + x[j][2] = x[i][2]; + v[j][0] = v[i][0]; + v[j][1] = v[i][1]; + v[j][2] = v[i][2]; + + q[j] = q[i]; + mu[j][0] = mu[i][0]; + mu[j][1] = mu[i][1]; + mu[j][2] = mu[i][2]; + omega[j][0] = omega[i][0]; + omega[j][1] = omega[i][1]; + omega[j][2] = omega[i][2]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j); +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecDipole::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = mu[j][0]; + buf[m++] = mu[j][1]; + buf[m++] = mu[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = mu[j][0]; + buf[m++] = mu[j][1]; + buf[m++] = mu[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecDipole::pack_comm_one(int i, double *buf) +{ + buf[0] = mu[i][0]; + buf[1] = mu[i][1]; + buf[2] = mu[i][2]; + return 3; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecDipole::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + mu[i][0] = buf[m++]; + mu[i][1] = buf[m++]; + mu[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecDipole::unpack_comm_one(int i, double *buf) +{ + mu[i][0] = buf[0]; + mu[i][1] = buf[1]; + mu[i][2] = buf[2]; + return 3; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecDipole::pack_reverse(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = f[i][0]; + buf[m++] = f[i][1]; + buf[m++] = f[i][2]; + buf[m++] = torque[i][0]; + buf[m++] = torque[i][1]; + buf[m++] = torque[i][2]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecDipole::pack_reverse_one(int i, double *buf) +{ + buf[0] = torque[i][0]; + buf[1] = torque[i][1]; + buf[2] = torque[i][2]; + return 3; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecDipole::unpack_reverse(int n, int *list, double *buf) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + f[j][0] += buf[m++]; + f[j][1] += buf[m++]; + f[j][2] += buf[m++]; + torque[j][0] += buf[m++]; + torque[j][1] += buf[m++]; + torque[j][2] += buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecDipole::unpack_reverse_one(int i, double *buf) +{ + torque[i][0] = buf[0]; + torque[i][1] = buf[1]; + torque[i][2] = buf[2]; + return 3; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecDipole::pack_border(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = mu[j][0]; + buf[m++] = mu[j][1]; + buf[m++] = mu[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = mu[j][0]; + buf[m++] = mu[j][1]; + buf[m++] = mu[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecDipole::pack_border_one(int i, double *buf) +{ + buf[0] = q[i]; + buf[1] = mu[i][0]; + buf[2] = mu[i][1]; + buf[3] = mu[i][2]; + return 4; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecDipole::unpack_border(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + q[i] = buf[m++]; + mu[i][0] = buf[m++]; + mu[i][1] = buf[m++]; + mu[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecDipole::unpack_border_one(int i, double *buf) +{ + q[i] = buf[0]; + mu[i][0] = buf[1]; + mu[i][1] = buf[2]; + mu[i][2] = buf[3]; + return 4; +} + +/* ---------------------------------------------------------------------- + pack all atom quantities for shipping to another proc + xyz must be 1st 3 values, so that comm::exchange can test on them +------------------------------------------------------------------------- */ + +int AtomVecDipole::pack_exchange(int i, double *buf) +{ + int m = 1; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + buf[m++] = v[i][0]; + buf[m++] = v[i][1]; + buf[m++] = v[i][2]; + buf[m++] = tag[i]; + buf[m++] = type[i]; + buf[m++] = mask[i]; + buf[m++] = image[i]; + + buf[m++] = q[i]; + buf[m++] = mu[i][0]; + buf[m++] = mu[i][1]; + buf[m++] = mu[i][2]; + buf[m++] = omega[i][0]; + buf[m++] = omega[i][1]; + buf[m++] = omega[i][2]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]); + + buf[0] = m; + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecDipole::unpack_exchange(double *buf) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + int m = 1; + x[nlocal][0] = buf[m++]; + x[nlocal][1] = buf[m++]; + x[nlocal][2] = buf[m++]; + v[nlocal][0] = buf[m++]; + v[nlocal][1] = buf[m++]; + v[nlocal][2] = buf[m++]; + tag[nlocal] = static_cast (buf[m++]); + type[nlocal] = static_cast (buf[m++]); + mask[nlocal] = static_cast (buf[m++]); + image[nlocal] = static_cast (buf[m++]); + + q[nlocal] = buf[m++]; + mu[nlocal][0] = buf[m++]; + mu[nlocal][1] = buf[m++]; + mu[nlocal][2] = buf[m++]; + omega[nlocal][0] = buf[m++]; + omega[nlocal][1] = buf[m++]; + omega[nlocal][2] = buf[m++]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + m += modify->fix[atom->extra_grow[iextra]]-> + unpack_exchange(nlocal,&buf[m]); + + atom->nlocal++; + return m; +} + +/* ---------------------------------------------------------------------- + size of restart data for all atoms owned by this proc + include extra data stored by fixes +------------------------------------------------------------------------- */ + +int AtomVecDipole::size_restart() +{ + int i; + + int nlocal = atom->nlocal; + int n = 18 * nlocal; + + if (atom->nextra_restart) + for (int iextra = 0; iextra < atom->nextra_restart; iextra++) + for (i = 0; i < nlocal; i++) + n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); + + return n; +} + +/* ---------------------------------------------------------------------- + pack atom I's data for restart file including extra quantities + xyz must be 1st 3 values, so that read_restart can test on them + molecular types may be negative, but write as positive +------------------------------------------------------------------------- */ + +int AtomVecDipole::pack_restart(int i, double *buf) +{ + int m = 1; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + buf[m++] = tag[i]; + buf[m++] = type[i]; + buf[m++] = mask[i]; + buf[m++] = image[i]; + buf[m++] = v[i][0]; + buf[m++] = v[i][1]; + buf[m++] = v[i][2]; + + buf[m++] = q[i]; + buf[m++] = mu[i][0]; + buf[m++] = mu[i][1]; + buf[m++] = mu[i][2]; + buf[m++] = omega[i][0]; + buf[m++] = omega[i][1]; + buf[m++] = omega[i][2]; + + if (atom->nextra_restart) + for (int iextra = 0; iextra < atom->nextra_restart; iextra++) + m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]); + + buf[0] = m; + return m; +} + +/* ---------------------------------------------------------------------- + unpack data for one atom from restart file including extra quantities +------------------------------------------------------------------------- */ + +int AtomVecDipole::unpack_restart(double *buf) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) { + grow(0); + if (atom->nextra_store) + atom->extra = memory->grow_2d_double_array(atom->extra,nmax, + atom->nextra_store, + "atom:extra"); + } + + int m = 1; + x[nlocal][0] = buf[m++]; + x[nlocal][1] = buf[m++]; + x[nlocal][2] = buf[m++]; + tag[nlocal] = static_cast (buf[m++]); + type[nlocal] = static_cast (buf[m++]); + mask[nlocal] = static_cast (buf[m++]); + image[nlocal] = static_cast (buf[m++]); + v[nlocal][0] = buf[m++]; + v[nlocal][1] = buf[m++]; + v[nlocal][2] = buf[m++]; + + q[nlocal] = buf[m++]; + mu[nlocal][0] = buf[m++]; + mu[nlocal][1] = buf[m++]; + mu[nlocal][2] = buf[m++]; + omega[nlocal][0] = buf[m++]; + omega[nlocal][1] = buf[m++]; + omega[nlocal][2] = buf[m++]; + + double **extra = atom->extra; + if (atom->nextra_store) { + int size = static_cast (buf[0]) - m; + for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++]; + } + + atom->nlocal++; + return m; +} + +/* ---------------------------------------------------------------------- + create one atom of itype at coord + set other values to defaults +------------------------------------------------------------------------- */ + +void AtomVecDipole::create_atom(int itype, double *coord) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + tag[nlocal] = 0; + type[nlocal] = itype; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + mask[nlocal] = 1; + image[nlocal] = (512 << 20) | (512 << 10) | 512; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + + q[nlocal] = 0.0; + mu[nlocal][0] = 0.0; + mu[nlocal][1] = 0.0; + mu[nlocal][2] = 0.0; + omega[nlocal][0] = 0.0; + omega[nlocal][1] = 0.0; + omega[nlocal][2] = 0.0; + + atom->nlocal++; +} + +/* ---------------------------------------------------------------------- + unpack one line from Atoms section of data file + initialize other atom quantities +------------------------------------------------------------------------- */ + +void AtomVecDipole::data_atom(double *coord, int imagetmp, char **values) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + tag[nlocal] = atoi(values[0]); + if (tag[nlocal] <= 0) + error->one("Invalid atom ID in Atoms section of data file"); + + type[nlocal] = atoi(values[1]); + if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) + error->one("Invalid atom type in Atoms section of data file"); + + q[nlocal] = atof(values[2]); + + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + + mu[nlocal][0] = atof(values[6]); + mu[nlocal][1] = atof(values[7]); + mu[nlocal][2] = atof(values[8]); + + image[nlocal] = imagetmp; + + mask[nlocal] = 1; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + omega[nlocal][0] = 0.0; + omega[nlocal][1] = 0.0; + omega[nlocal][2] = 0.0; + + atom->nlocal++; +} + +/* ---------------------------------------------------------------------- + unpack hybrid quantities from one line in Atoms section of data file + initialize other atom quantities for this sub-style +------------------------------------------------------------------------- */ + +int AtomVecDipole::data_atom_hybrid(int nlocal, char **values) +{ + q[nlocal] = atof(values[0]); + mu[nlocal][0] = atof(values[1]); + mu[nlocal][1] = atof(values[2]); + mu[nlocal][2] = atof(values[3]); + + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + omega[nlocal][0] = 0.0; + omega[nlocal][1] = 0.0; + omega[nlocal][2] = 0.0; + + return 4; +} + +/* ---------------------------------------------------------------------- + unpack one line from Velocities section of data file +------------------------------------------------------------------------- */ + +void AtomVecDipole::data_vel(int m, char **values) +{ + v[m][0] = atof(values[0]); + v[m][1] = atof(values[1]); + v[m][2] = atof(values[2]); + omega[m][0] = atof(values[3]); + omega[m][1] = atof(values[4]); + omega[m][2] = atof(values[5]); +} + +/* ---------------------------------------------------------------------- + unpack hybrid quantities from one line in Velocities section of data file +------------------------------------------------------------------------- */ + +int AtomVecDipole::data_vel_hybrid(int m, char **values) +{ + omega[m][0] = atof(values[0]); + omega[m][1] = atof(values[1]); + omega[m][2] = atof(values[2]); + return 3; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory +------------------------------------------------------------------------- */ + +int AtomVecDipole::memory_usage() +{ + int bytes = 0; + + if (atom->memcheck("tag")) bytes += nmax * sizeof(int); + if (atom->memcheck("type")) bytes += nmax * sizeof(int); + if (atom->memcheck("mask")) bytes += nmax * sizeof(int); + if (atom->memcheck("image")) bytes += nmax * sizeof(int); + if (atom->memcheck("x")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("v")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("f")) bytes += nmax*3 * sizeof(double); + + if (atom->memcheck("q")) bytes += nmax * sizeof(double); + if (atom->memcheck("mu")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("omega")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("torque")) bytes += nmax*3 * sizeof(double); + + return bytes; +} diff --git a/src/DIPOLE/atom_vec_dipole.h b/src/DIPOLE/atom_vec_dipole.h new file mode 100644 index 0000000000..516dbc9223 --- /dev/null +++ b/src/DIPOLE/atom_vec_dipole.h @@ -0,0 +1,58 @@ +/* ---------------------------------------------------------------------- + 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 ATOM_VEC_DIPOLE_H +#define ATOM_VEC_DIPOLE_H + +#include "atom_vec.h" + +namespace LAMMPS_NS { + +class AtomVecDipole : public AtomVec { + public: + AtomVecDipole(class LAMMPS *, int, char **); + void grow(int); + void copy(int, int); + int pack_comm(int, int *, double *, int, int *); + int pack_comm_one(int, double *); + void unpack_comm(int, int, double *); + int unpack_comm_one(int, double *); + int pack_reverse(int, int, double *); + int pack_reverse_one(int, double *); + void unpack_reverse(int, int *, double *); + int unpack_reverse_one(int, double *); + int pack_border(int, int *, double *, int, int *); + int pack_border_one(int, double *); + void unpack_border(int, int, double *); + int unpack_border_one(int, double *); + int pack_exchange(int, double *); + int unpack_exchange(double *); + int size_restart(); + int pack_restart(int, double *); + int unpack_restart(double *); + void create_atom(int, double *); + void data_atom(double *, int, char **); + int data_atom_hybrid(int, char **); + void data_vel(int, char **); + int data_vel_hybrid(int, char **); + int memory_usage(); + + private: + int *tag,*type,*mask,*image; + double **x,**v,**f; + double *q,**mu,**omega,**torque; +}; + +} + +#endif diff --git a/src/DIPOLE/compute_temp_dipole.cpp b/src/DIPOLE/compute_temp_dipole.cpp new file mode 100644 index 0000000000..abb5d54bd1 --- /dev/null +++ b/src/DIPOLE/compute_temp_dipole.cpp @@ -0,0 +1,144 @@ +/* ---------------------------------------------------------------------- + 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 "group.h" +#include "modify.h" +#include "fix.h" +#include "error.h" + +using namespace LAMMPS_NS; + +// moment of inertia for a sphere + +#define INERTIA 0.4 + +/* ---------------------------------------------------------------------- */ + +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 == NULL || atom->shape == NULL) + error->all("Compute temp/dipole requires atom attributes omega, shape"); + + scalar_flag = vector_flag = 1; + size_vector = 6; + extensive = 0; + 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 * force->dimension * natoms; + dof -= extra_dof + fix_dof; + if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); + else tfactor = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +double ComputeTempDipole::compute_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; + + 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; +} diff --git a/src/DIPOLE/compute_temp_dipole.h b/src/DIPOLE/compute_temp_dipole.h new file mode 100644 index 0000000000..74fafd84d9 --- /dev/null +++ b/src/DIPOLE/compute_temp_dipole.h @@ -0,0 +1,39 @@ +/* ---------------------------------------------------------------------- + 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 diff --git a/src/DIPOLE/fix_nve_dipole.cpp b/src/DIPOLE/fix_nve_dipole.cpp new file mode 100644 index 0000000000..d056d324b6 --- /dev/null +++ b/src/DIPOLE/fix_nve_dipole.cpp @@ -0,0 +1,189 @@ +/* ---------------------------------------------------------------------- + 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 "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 == NULL || atom->omega == NULL || + atom->torque == NULL || atom->shape == NULL) + 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() +{ + 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; + + double g[3],h[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; + + // 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 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(); + else final_integrate(); +} + +/* ---------------------------------------------------------------------- */ + +void FixNVEDipole::final_integrate_respa(int ilevel) +{ + dtf = 0.5 * step_respa[ilevel] * force->ftm2v; + final_integrate(); +} diff --git a/src/DIPOLE/fix_nve_dipole.h b/src/DIPOLE/fix_nve_dipole.h new file mode 100644 index 0000000000..b19a0ff97b --- /dev/null +++ b/src/DIPOLE/fix_nve_dipole.h @@ -0,0 +1,40 @@ +/* ---------------------------------------------------------------------- + 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(); + void final_integrate(); + void initial_integrate_respa(int,int); + void final_integrate_respa(int); + + private: + double dtv,dtf; + double *step_respa; + double *inertia; +}; + +} + +#endif diff --git a/src/DIPOLE/pair_dipole_cut.cpp b/src/DIPOLE/pair_dipole_cut.cpp new file mode 100644 index 0000000000..28b470f24a --- /dev/null +++ b/src/DIPOLE/pair_dipole_cut.cpp @@ -0,0 +1,625 @@ +/* ---------------------------------------------------------------------- + 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 "stdlib.h" +#include "pair_dipole_cut.h" +#include "atom.h" +#include "neighbor.h" +#include "comm.h" +#include "force.h" +#include "update.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairDipoleCut::PairDipoleCut(LAMMPS *lmp) : Pair(lmp) {} + +/* ---------------------------------------------------------------------- */ + +PairDipoleCut::~PairDipoleCut() +{ + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + + memory->destroy_2d_double_array(cut_lj); + memory->destroy_2d_double_array(cut_ljsq); + memory->destroy_2d_double_array(cut_coul); + memory->destroy_2d_double_array(cut_coulsq); + memory->destroy_2d_double_array(epsilon); + memory->destroy_2d_double_array(sigma); + memory->destroy_2d_double_array(lj1); + memory->destroy_2d_double_array(lj2); + memory->destroy_2d_double_array(lj3); + memory->destroy_2d_double_array(lj4); + memory->destroy_2d_double_array(offset); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairDipoleCut::compute(int eflag, int vflag) +{ + int i,j,k,numneigh,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz; + double rsq,rinv,r2inv,r6inv,r3inv,r5inv,r7inv; + double forcecoulx,forcecouly,forcecoulz,fforce,crossx,crossy,crossz; + double tixcoul,tiycoul,tizcoul,tjxcoul,tjycoul,tjzcoul; + double fq,fx,fy,fz; + double pdotp,pidotr,pjdotr,pre1,pre2,pre3,pre4; + double forcelj,factor_coul,factor_lj; + double factor,phicoul,philj; + int *neighs; + double **f; + + eng_vdwl = eng_coul = 0.0; + if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; + + if (vflag == 2) f = update->f_pair; + else f = atom->f; + double **x = atom->x; + double *q = atom->q; + double **mu = atom->mu; + double **torque = atom->torque; + int *type = atom->type; + double *dipole = atom->dipole; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + // loop over neighbors of my atoms + + for (i = 0; i < nlocal; i++) { + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + neighs = neighbor->firstneigh[i]; + numneigh = neighbor->numneigh[i]; + + for (k = 0; k < numneigh; k++) { + j = neighs[k]; + + if (j < nall) factor_coul = factor_lj = 1.0; + else { + factor_coul = special_coul[j/nall]; + factor_lj = special_lj[j/nall]; + j = j % nall; + } + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r2inv = 1.0/rsq; + rinv = sqrt(r2inv); + + // atom can have both a charge and dipole + // i,j = charge-charge, dipole-dipole, dipole-charge, or charge-dipole + + if (rsq < cut_coulsq[itype][jtype]) { + + forcecoulx = forcecouly = forcecoulz = 0.0; + tixcoul = tiycoul = tizcoul = 0.0; + tjxcoul = tjycoul = tjzcoul = 0.0; + + if (qtmp != 0.0 && q[j] != 0.0) { + r3inv = r2inv*rinv; + pre1 = qtmp*q[j]*r3inv; + + forcecoulx += pre1*delx; + forcecouly += pre1*dely; + forcecoulz += pre1*delz; + } + + if (dipole[itype] > 0.0 && dipole[jtype] > 0.0) { + r3inv = r2inv*rinv; + r5inv = r3inv*r2inv; + r7inv = r5inv*r2inv; + + pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2]; + pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; + pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; + + pre1 = 3.0*r5inv*pdotp - 15.0*r7inv*pidotr*pjdotr; + pre2 = 3.0*r5inv*pjdotr; + pre3 = 3.0*r5inv*pidotr; + pre4 = -1.0*r3inv; + + forcecoulx += pre1*delx + pre2*mu[i][0] + pre3*mu[j][0]; + forcecouly += pre1*dely + pre2*mu[i][1] + pre3*mu[j][1]; + forcecoulz += pre1*delz + pre2*mu[i][2] + pre3*mu[j][2]; + + crossx = pre4 * (mu[i][1]*mu[j][2] - mu[i][2]*mu[j][1]); + crossy = pre4 * (mu[i][2]*mu[j][0] - mu[i][0]*mu[j][2]); + crossz = pre4 * (mu[i][0]*mu[j][1] - mu[i][1]*mu[j][0]); + + tixcoul += crossx + pre2 * (mu[i][1]*delz - mu[i][2]*dely); + tiycoul += crossy + pre2 * (mu[i][2]*delx - mu[i][0]*delz); + tizcoul += crossz + pre2 * (mu[i][0]*dely - mu[i][1]*delx); + tjxcoul += -crossx + pre3 * (mu[j][1]*delz - mu[j][2]*dely); + tjycoul += -crossy + pre3 * (mu[j][2]*delx - mu[j][0]*delz); + tjzcoul += -crossz + pre3 * (mu[j][0]*dely - mu[j][1]*delx); + } + + if (dipole[itype] > 0.0 && q[j] != 0.0) { + r3inv = r2inv*rinv; + r5inv = r3inv*r2inv; + pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; + pre1 = 3.0*q[j]*r5inv * pidotr; + pre2 = q[j]*r3inv; + + forcecoulx += pre2*mu[i][0] - pre1*delx; + forcecouly += pre2*mu[i][1] - pre1*dely; + forcecoulz += pre2*mu[i][2] - pre1*delz; + tixcoul += pre2 * (mu[i][1]*delz - mu[i][2]*dely); + tiycoul += pre2 * (mu[i][2]*delx - mu[i][0]*delz); + tizcoul += pre2 * (mu[i][0]*dely - mu[i][1]*delx); + } + + if (dipole[jtype] > 0.0 && qtmp != 0.0) { + r3inv = r2inv*rinv; + r5inv = r3inv*r2inv; + pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; + pre1 = 3.0*qtmp*r5inv * pjdotr; + pre2 = qtmp*r3inv; + + forcecoulx += pre1*delx - pre2*mu[j][0]; + forcecouly += pre1*dely - pre2*mu[j][1]; + forcecoulz += pre1*delz - pre2*mu[j][2]; + tjxcoul += -pre2 * (mu[j][1]*delz - mu[j][2]*dely); + tjycoul += -pre2 * (mu[j][2]*delx - mu[j][0]*delz); + tjzcoul += -pre2 * (mu[j][0]*dely - mu[j][1]*delx); + } + } + + // LJ interaction + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + fforce = factor_lj * forcelj*r2inv; + } else fforce = 0.0; + + // total force + + fq = factor_coul*qqrd2e; + fx = fq*forcecoulx + delx*fforce; + fy = fq*forcecouly + dely*fforce; + fz = fq*forcecoulz + delz*fforce; + + // force & torque accumulation + + f[i][0] += fx; + f[i][1] += fy; + f[i][2] += fz; + torque[i][0] += fq*tixcoul; + torque[i][1] += fq*tiycoul; + torque[i][2] += fq*tizcoul; + + if (newton_pair || j < nlocal) { + f[j][0] -= fx; + f[j][1] -= fy; + f[j][2] -= fz; + torque[j][0] += fq*tjxcoul; + torque[j][1] += fq*tjycoul; + torque[j][2] += fq*tjzcoul; + } + + if (eflag) { + if (newton_pair || j < nlocal) factor = 1.0; + else factor = 0.5; + + if (rsq < cut_coulsq[itype][jtype]) { + phicoul = qtmp*q[j]*rinv; + if (dipole[itype] > 0.0 && dipole[jtype] > 0.0) + phicoul += r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr; + if (dipole[itype] > 0.0 && q[j] != 0.0) + phicoul += -q[j]*r3inv*pidotr; + if (dipole[jtype] > 0.0 && qtmp != 0.0) + phicoul += qtmp*r3inv*pjdotr; + eng_coul += factor*factor_coul*qqrd2e*phicoul; + } + + if (rsq < cut_ljsq[itype][jtype]) { + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + eng_vdwl += factor*factor_lj*philj; + } + } + + if (vflag == 1) { + if (newton_pair == 0 && j >= nlocal) { + fx *= 0.5; fy *= 0.5; fz *= 0.5; + } + virial[0] += delx*fx; + virial[1] += dely*fy; + virial[2] += delz*fz; + virial[3] += delx*fy; + virial[4] += delx*fz; + virial[5] += dely*fz; + } + } + } + } + if (vflag == 2) virial_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairDipoleCut::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + + cut_lj = memory->create_2d_double_array(n+1,n+1,"pair:cut_lj"); + cut_ljsq = memory->create_2d_double_array(n+1,n+1,"pair:cut_ljsq"); + cut_coul = memory->create_2d_double_array(n+1,n+1,"pair:cut_coul"); + cut_coulsq = memory->create_2d_double_array(n+1,n+1,"pair:cut_coulsq"); + epsilon = memory->create_2d_double_array(n+1,n+1,"pair:epsilon"); + sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma"); + lj1 = memory->create_2d_double_array(n+1,n+1,"pair:lj1"); + lj2 = memory->create_2d_double_array(n+1,n+1,"pair:lj2"); + lj3 = memory->create_2d_double_array(n+1,n+1,"pair:lj3"); + lj4 = memory->create_2d_double_array(n+1,n+1,"pair:lj4"); + offset = memory->create_2d_double_array(n+1,n+1,"pair:offset"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairDipoleCut::settings(int narg, char **arg) +{ + if (narg < 1 || narg > 2) + error->all("Incorrect args in pair_style command"); + + cut_lj_global = atof(arg[0]); + if (narg == 1) cut_coul_global = cut_lj_global; + else cut_coul_global = atof(arg[1]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) { + cut_lj[i][j] = cut_lj_global; + cut_coul[i][j] = cut_coul_global; + } + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairDipoleCut::coeff(int narg, char **arg) +{ + if (narg < 4 || narg > 6) + error->all("Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double epsilon_one = atof(arg[2]); + double sigma_one = atof(arg[3]); + + double cut_lj_one = cut_lj_global; + double cut_coul_one = cut_coul_global; + if (narg >= 5) cut_coul_one = cut_lj_one = atof(arg[4]); + if (narg == 6) cut_coul_one = atof(arg[5]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + epsilon[i][j] = epsilon_one; + sigma[i][j] = sigma_one; + cut_lj[i][j] = cut_lj_one; + cut_coul[i][j] = cut_coul_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairDipoleCut::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], + sigma[i][i],sigma[j][j]); + sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); + cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); + cut_coul[i][j] = mix_distance(cut_coul[i][i],cut_coul[j][j]); + } + + double cut = MAX(cut_lj[i][j],cut_coul[i][j]); + cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; + cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j]; + + lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + + if (offset_flag) { + double ratio = sigma[i][j] / cut_lj[i][j]; + offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); + } else offset[i][j] = 0.0; + + cut_ljsq[j][i] = cut_ljsq[i][j]; + cut_coulsq[j][i] = cut_coulsq[i][j]; + lj1[j][i] = lj1[i][j]; + lj2[j][i] = lj2[i][j]; + lj3[j][i] = lj3[i][j]; + lj4[j][i] = lj4[i][j]; + offset[j][i] = offset[i][j]; + + return cut; +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairDipoleCut::init_style() +{ + if (atom->q == NULL || atom->mu == NULL || + atom->torque == NULL || atom->dipole == NULL) + error->all("Pair dipole/cut requires atom attributes " + "q, mu, torque, dipole"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairDipoleCut::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&epsilon[i][j],sizeof(double),1,fp); + fwrite(&sigma[i][j],sizeof(double),1,fp); + fwrite(&cut_lj[i][j],sizeof(double),1,fp); + fwrite(&cut_coul[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairDipoleCut::read_restart(FILE *fp) +{ + read_restart_settings(fp); + + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&epsilon[i][j],sizeof(double),1,fp); + fread(&sigma[i][j],sizeof(double),1,fp); + fread(&cut_lj[i][j],sizeof(double),1,fp); + fread(&cut_coul[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairDipoleCut::write_restart_settings(FILE *fp) +{ + fwrite(&cut_lj_global,sizeof(double),1,fp); + fwrite(&cut_coul_global,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairDipoleCut::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_lj_global,sizeof(double),1,fp); + fread(&cut_coul_global,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +void PairDipoleCut::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, int eflag, + One &one) +{ + double rinv,r2inv,r6inv,r3inv,r5inv,r7inv; + double forcecoulx,forcecouly,forcecoulz,fforce,crossx,crossy,crossz; + double tixcoul,tiycoul,tizcoul,tjxcoul,tjycoul,tjzcoul; + double pdotp,pdotr,pidotr,pjdotr,pre1,pre2,pre3,pre4; + double fq,forcelj,phicoul,philj; + + double delx = atom->x[i][0] - atom->x[j][0]; + double dely = atom->x[i][1] - atom->x[j][1]; + double delz = atom->x[i][2] - atom->x[j][2]; + + r2inv = 1.0/rsq; + rinv = sqrt(r2inv); + + forcecoulx = forcecouly = forcecoulz = 0.0; + tixcoul = tiycoul = tizcoul = 0.0; + tjxcoul = tjycoul = tjzcoul = 0.0; + + if (rsq < cut_coulsq[itype][jtype]) { + double **mu = atom->mu; + + if (atom->q[i] != 0.0 && atom->q[j] != 0.0) { + r3inv = r2inv*rinv; + pre1 = atom->q[i]*atom->q[j]*r3inv; + + forcecoulx += pre1*delx; + forcecouly += pre1*dely; + forcecoulz += pre1*delz; + } + + if (atom->dipole[itype] > 0.0 && atom->dipole[jtype] > 0.0) { + r3inv = r2inv*rinv; + r5inv = r3inv*r2inv; + r7inv = r5inv*r2inv; + + pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2]; + pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; + pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; + + pre1 = 3.0*r5inv*pdotp - 15.0*r7inv*pidotr*pjdotr; + pre2 = 3.0*r5inv*pjdotr; + pre3 = 3.0*r5inv*pidotr; + pre4 = -1.0*r3inv; + + forcecoulx += pre1*delx + pre2*mu[i][0] + pre3*mu[j][0]; + forcecouly += pre1*dely + pre2*mu[i][1] + pre3*mu[j][1]; + forcecoulz += pre1*delz + pre2*mu[i][2] + pre3*mu[j][2]; + + crossx = pre4 * (mu[i][1]*mu[j][2] - mu[i][2]*mu[j][1]); + crossy = pre4 * (mu[i][2]*mu[j][0] - mu[i][0]*mu[j][2]); + crossz = pre4 * (mu[i][0]*mu[j][1] - mu[i][1]*mu[j][0]); + + tixcoul += crossx + pre2 * (mu[i][1]*delz - mu[i][2]*dely); + tiycoul += crossy + pre2 * (mu[i][2]*delx - mu[i][0]*delz); + tizcoul += crossz + pre2 * (mu[i][0]*dely - mu[i][1]*delx); + tjxcoul += -crossx + pre3 * (mu[j][1]*delz - mu[j][2]*dely); + tjycoul += -crossy + pre3 * (mu[j][2]*delx - mu[j][0]*delz); + tjzcoul += -crossz + pre3 * (mu[j][0]*dely - mu[j][1]*delx); + + } else if (atom->dipole[itype] > 0.0 && atom->q[j] != 0.0) { + r3inv = r2inv*rinv; + r5inv = r3inv*r2inv; + pdotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; + pre1 = 3.0*atom->q[j]*r5inv * pdotr; + pre2 = atom->q[j]*r3inv; + + forcecoulx += pre2*mu[i][0] - pre1*delx; + forcecouly += pre2*mu[i][1] - pre1*dely; + forcecoulz += pre2*mu[i][2] - pre1*delz; + tixcoul += pre2 * (mu[i][1]*delz - mu[i][2]*dely); + tiycoul += pre2 * (mu[i][2]*delx - mu[i][0]*delz); + tizcoul += pre2 * (mu[i][0]*dely - mu[i][1]*delx); + + } else if (atom->dipole[jtype] > 0.0 && atom->q[i] != 0.0) { + r3inv = r2inv*rinv; + r5inv = r3inv*r2inv; + pdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; + pre1 = 3.0*atom->q[i]*r5inv * pdotr; + pre2 = atom->q[i]*r3inv; + + forcecoulx += pre1*delx - pre2*mu[j][0]; + forcecouly += pre1*dely - pre2*mu[j][1]; + forcecoulz += pre1*delz - pre2*mu[j][2]; + tjxcoul += -pre2 * (mu[j][1]*delz - mu[j][2]*dely); + tjycoul += -pre2 * (mu[j][2]*delx - mu[j][0]*delz); + tjzcoul += -pre2 * (mu[j][0]*dely - mu[j][1]*delx); + } + } + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + fforce = factor_lj * forcelj*r2inv; + } else fforce = 0.0; + + fq = factor_coul*force->qqrd2e; + one.fx = fq*forcecoulx + delx*fforce; + one.fy = fq*forcecouly + dely*fforce; + one.fz = fq*forcecoulz + delz*fforce; + one.tix = fq*tixcoul; + one.tiy = fq*tiycoul; + one.tiz = fq*tizcoul; + one.tjx = fq*tjxcoul; + one.tjy = fq*tjycoul; + one.tjz = fq*tjzcoul; + + if (eflag) { + if (rsq < cut_coulsq[itype][jtype]) { + phicoul = atom->q[i]*atom->q[j]*rinv; + if (atom->dipole[itype] > 0.0 && atom->dipole[jtype] > 0.0) + phicoul += r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr; + else if (atom->dipole[itype] > 0.0 && atom->q[j] != 0.0) + phicoul += -pre2*pdotr; + else if (atom->dipole[jtype] > 0.0 && atom->q[i] != 0.0) + phicoul += pre2*pdotr; + one.eng_coul = factor_coul*force->qqrd2e*phicoul; + } else one.eng_coul = 0.0; + if (rsq < cut_ljsq[itype][jtype]) { + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + one.eng_vdwl = factor_lj*philj; + } else one.eng_vdwl = 0.0; + } +} diff --git a/src/DIPOLE/pair_dipole_cut.h b/src/DIPOLE/pair_dipole_cut.h new file mode 100644 index 0000000000..e3a5367405 --- /dev/null +++ b/src/DIPOLE/pair_dipole_cut.h @@ -0,0 +1,48 @@ +/* ---------------------------------------------------------------------- + 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 PAIR_DIPOLE_CUT_H +#define PAIR_DIPOLE_CUT_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairDipoleCut : public Pair { + public: + PairDipoleCut(class LAMMPS *); + ~PairDipoleCut(); + void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void init_style(); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + void single(int, int, int, int, double, double, double, int, One &); + + private: + double cut_lj_global,cut_coul_global; + double **cut_lj,**cut_ljsq; + double **cut_coul,**cut_coulsq; + double **epsilon,**sigma; + double **lj1,**lj2,**lj3,**lj4,**offset; + + void allocate(); +}; + +} + +#endif diff --git a/src/DIPOLE/style_dipole.h b/src/DIPOLE/style_dipole.h new file mode 100644 index 0000000000..7f989773f4 --- /dev/null +++ b/src/DIPOLE/style_dipole.h @@ -0,0 +1,44 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef AtomInclude +#include "atom_vec_dipole.h" +#endif + +#ifdef AtomClass +AtomStyle(dipole,AtomVecDipole) +#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 +#include "pair_dipole_cut.h" +#endif + +#ifdef PairClass +PairStyle(dipole/cut,PairDipoleCut) +#endif diff --git a/src/MANYBODY/Install.csh b/src/MANYBODY/Install.csh index 5cfed7a791..3ac57d2ba8 100644 --- a/src/MANYBODY/Install.csh +++ b/src/MANYBODY/Install.csh @@ -4,14 +4,12 @@ if ($1 == 1) then cp style_manybody.h .. - cp pair_airebo.cpp .. cp pair_eam.cpp .. cp pair_eam_alloy.cpp .. cp pair_eam_fs.cpp .. cp pair_sw.cpp .. cp pair_tersoff.cpp .. - cp pair_airebo.h .. cp pair_eam.h .. cp pair_eam_alloy.h .. cp pair_eam_fs.h .. @@ -23,14 +21,12 @@ else if ($1 == 0) then rm ../style_manybody.h touch ../style_manybody.h - rm ../pair_airebo.cpp rm ../pair_eam.cpp rm ../pair_eam_alloy.cpp rm ../pair_eam_fs.cpp rm ../pair_sw.cpp rm ../pair_tersoff.cpp - rm ../pair_airebo.h rm ../pair_eam.h rm ../pair_eam_alloy.h rm ../pair_eam_fs.h diff --git a/src/MANYBODY/style_manybody.h b/src/MANYBODY/style_manybody.h index c564993e3e..449895a530 100644 --- a/src/MANYBODY/style_manybody.h +++ b/src/MANYBODY/style_manybody.h @@ -12,7 +12,6 @@ ------------------------------------------------------------------------- */ #ifdef PairInclude -#include "pair_airebo.h" #include "pair_eam.h" #include "pair_eam_alloy.h" #include "pair_eam_fs.h" @@ -21,7 +20,6 @@ #endif #ifdef PairClass -PairStyle(airebo,PairAIREBO) PairStyle(eam,PairEAM) PairStyle(eam/alloy,PairEAMAlloy) PairStyle(eam/fs,PairEAMFS) diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp new file mode 100644 index 0000000000..35831ec4f0 --- /dev/null +++ b/src/compute_temp_deform.cpp @@ -0,0 +1,193 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Pieter in 't Veld (SNL) +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "string.h" +#include "compute_temp_deform.h" +#include "domain.h" +#include "atom.h" +#include "force.h" +#include "modify.h" +#include "fix.h" +#include "fix_deform.h" +#include "group.h" +#include "comm.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp + +/* ---------------------------------------------------------------------- */ + +ComputeTempDeform::ComputeTempDeform(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 3) error->all("Illegal compute temp/deform command"); + + scalar_flag = vector_flag = 1; + size_vector = 6; + extensive = 0; + tempflag = 1; + + vector = new double[6]; +} + +/* ---------------------------------------------------------------------- */ + +ComputeTempDeform::~ComputeTempDeform() +{ + delete [] vector; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempDeform::init() +{ + int i; + + fix_dof = 0; + for (i = 0; i < modify->nfix; i++) + fix_dof += modify->fix[i]->dof(igroup); + recount(); + + // check fix deform remap settings + + for (i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"deform") == 0) { + if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP && + comm->me == 0) + error->warning("Using fix nvt/sllod with inconsistent fix deform remapping"); + break; + } + if (i == modify->nfix && comm->me == 0) + error->warning("Using fix nvt/sllod with no fix deform defined"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempDeform::recount() +{ + double natoms = group->count(igroup); + dof = force->dimension * natoms; + dof -= extra_dof + fix_dof; + if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); + else tfactor = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +double ComputeTempDeform::compute_scalar() +{ + double lamda[3],vstream[3],vthermal[3]; + + double **x = atom->x; + double **v = atom->v; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + // lamda = 0-1 triclinic lamda coords + // vstream = streaming velocity = Hrate*lamda + Hratelo + // vthermal = thermal velocity = v - vstream + + double *h_rate = domain->h_rate; + double *h_ratelo = domain->h_ratelo; + + double t = 0.0; + + if (mass) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + vthermal[0] = v[i][0] - vstream[0]; + vthermal[1] = v[i][1] - vstream[1]; + vthermal[2] = v[i][2] - vstream[2]; + t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + + vthermal[2]*vthermal[2]) * mass[type[i]]; + } + } + else + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + vthermal[0] = v[i][0] - vstream[0]; + vthermal[1] = v[i][1] - vstream[1]; + vthermal[2] = v[i][2] - vstream[2]; + t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + + vthermal[2]*vthermal[2]) * rmass[i]; + } + + MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); + scalar *= tfactor; + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempDeform::compute_vector() +{ + double lamda[3],vstream[3],vthermal[3]; + + double **x = atom->x; + double **v = atom->v; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *h_rate = domain->h_rate; + double *h_ratelo = domain->h_ratelo; + + double massone,t[6]; + for (int i = 0; i < 6; i++) t[i] = 0.0; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + vthermal[0] = v[i][0] - vstream[0]; + vthermal[1] = v[i][1] - vstream[1]; + vthermal[2] = v[i][2] - vstream[2]; + + if (mass) massone = mass[type[i]]; + else massone = rmass[i]; + t[0] += massone * vthermal[0]*vthermal[0]; + t[1] += massone * vthermal[1]*vthermal[1]; + t[2] += massone * vthermal[2]*vthermal[2]; + t[3] += massone * vthermal[0]*vthermal[1]; + t[4] += massone * vthermal[0]*vthermal[2]; + t[5] += massone * vthermal[1]*vthermal[2]; + } + + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); + for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; +} diff --git a/src/compute_temp_deform.h b/src/compute_temp_deform.h new file mode 100644 index 0000000000..2eba8ed091 --- /dev/null +++ b/src/compute_temp_deform.h @@ -0,0 +1,38 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef COMPUTE_TEMP_DEFORM_H +#define COMPUTE_TEMP_DEFORM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeTempDeform : public Compute { + public: + ComputeTempDeform(class LAMMPS *, int, char **); + ~ComputeTempDeform(); + void init(); + double compute_scalar(); + void compute_vector(); + + private: + int fix_dof; + double tfactor; + + void recount(); +}; + +} + +#endif diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp new file mode 100644 index 0000000000..6bfe803465 --- /dev/null +++ b/src/fix_deform.cpp @@ -0,0 +1,674 @@ +/* ---------------------------------------------------------------------- + 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 author: Pieter in 't Veld (SNL) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "stdlib.h" +#include "math.h" +#include "fix_deform.h" +#include "atom.h" +#include "update.h" +#include "comm.h" +#include "domain.h" +#include "lattice.h" +#include "force.h" +#include "modify.h" +#include "kspace.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{NONE,FINAL,DELTA,SCALE,VEL,RATE,VOLUME}; +enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE}; + +// same as domain.cpp, fix_nvt_sllod.cpp, compute_temp_deform.cpp + +enum{NO_REMAP,X_REMAP,V_REMAP}; + +/* ---------------------------------------------------------------------- */ + +FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +{ + if (narg < 4) error->all("Illegal fix deform command"); + + box_change = 1; + + nevery = atoi(arg[3]); + if (nevery <= 0) error->all("Illegal fix deform command"); + + // set defaults + + set = new Set[6]; + set[0].style = set[1].style = set[2].style = + set[3].style = set[4].style = set[5].style = NONE; + + // parse arguments + + triclinic = domain->triclinic; + + int index; + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"x") == 0 || strcmp(arg[iarg],"y") == 0 || + strcmp(arg[iarg],"z") == 0) { + if (strcmp(arg[iarg],"x") == 0) index = 0; + else if (strcmp(arg[iarg],"y") == 0) index = 1; + else if (strcmp(arg[iarg],"z") == 0) index = 2; + + if (iarg+2 > narg) error->all("Illegal fix deform command"); + if (strcmp(arg[iarg+1],"final") == 0) { + if (iarg+4 > narg) error->all("Illegal fix deform command"); + set[index].style = FINAL; + set[index].flo = atof(arg[iarg+2]); + set[index].fhi = atof(arg[iarg+3]); + iarg += 4; + } else if (strcmp(arg[iarg+1],"delta") == 0) { + if (iarg+4 > narg) error->all("Illegal fix deform command"); + set[index].style = DELTA; + set[index].dlo = atof(arg[iarg+2]); + set[index].dhi = atof(arg[iarg+3]); + iarg += 4; + } else if (strcmp(arg[iarg+1],"scale") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deform command"); + set[index].style = SCALE; + set[index].scale = atof(arg[iarg+2]); + iarg += 3; + } else if (strcmp(arg[iarg+1],"vel") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deform command"); + set[index].style = VEL; + set[index].vel = atof(arg[iarg+2]); + iarg += 3; + } else if (strcmp(arg[iarg+1],"rate") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deform command"); + set[index].style = RATE; + set[index].rate = atof(arg[iarg+2]); + iarg += 3; + } else if (strcmp(arg[iarg+1],"volume") == 0) { + set[index].style = VOLUME; + iarg += 2; + } else error->all("Illegal fix deform command"); + + } else if (strcmp(arg[iarg],"xy") == 0 || strcmp(arg[iarg],"xz") == 0 || + strcmp(arg[iarg],"yz") == 0) { + if (triclinic == 0) + error->all("Fix deform tilt factors require triclinic box"); + if (strcmp(arg[iarg],"xy") == 0) index = 5; + else if (strcmp(arg[iarg],"xz") == 0) index = 4; + else if (strcmp(arg[iarg],"yz") == 0) index = 3; + if (iarg+2 > narg) error->all("Illegal fix deform command"); + if (strcmp(arg[iarg+1],"final") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deform command"); + set[index].style = FINAL; + set[index].ftilt = atof(arg[iarg+2]); + iarg += 3; + } else if (strcmp(arg[iarg+1],"delta") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deform command"); + set[index].style = DELTA; + set[index].dtilt = atof(arg[iarg+2]); + iarg += 3; + } else if (strcmp(arg[iarg+1],"vel") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deform command"); + set[index].style = VEL; + set[index].vel = atof(arg[iarg+2]); + iarg += 3; + } else if (strcmp(arg[iarg+1],"rate") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deform command"); + set[index].style = RATE; + set[index].rate = atof(arg[iarg+2]); + iarg += 3; + } else error->all("Illegal fix deform command"); + } else break; + } + + // read options from end of input line + + options(narg-iarg,&arg[iarg]); + + // check periodicity + + if ((set[0].style && domain->xperiodic == 0) || + (set[1].style && domain->yperiodic == 0) || + (set[2].style && domain->zperiodic == 0)) + error->all("Cannot fix deform on a non-periodic boundary"); + + if (set[3].style && (domain->yperiodic == 0 || domain->zperiodic == 0)) + error->all("Cannot fix deform on a non-periodic boundary"); + if (set[4].style && (domain->xperiodic == 0 || domain->zperiodic == 0)) + error->all("Cannot fix deform on a non-periodic boundary"); + if (set[5].style && (domain->xperiodic == 0 || domain->yperiodic == 0)) + error->all("Cannot fix deform on a non-periodic boundary"); + + // setup scaling + + if (scaleflag && domain->lattice == NULL) + error->all("Use of fix deform with undefined lattice"); + + if (scaleflag) { + xscale = domain->lattice->xlattice; + yscale = domain->lattice->ylattice; + zscale = domain->lattice->zlattice; + } + else xscale = yscale = zscale = 1.0; + + // apply scaling to input parameters with dist/vel units + // for 3,4,5 scale is in 1st dimension, e.g. x for xz + + double map[6]; + map[0] = xscale; map[1] = yscale; map[2] = zscale; + map[3] = yscale; map[4] = xscale; map[5] = xscale; + + for (int i = 0; i < 3; i++) { + if (set[i].style == FINAL) { + set[i].flo *= map[i]; + set[i].fhi *= map[i]; + } else if (set[i].style == DELTA) { + set[i].dlo *= map[i]; + set[i].dhi *= map[i]; + } else if (set[i].style == VEL) { + set[i].vel *= map[i]; + } + } + + for (int i = 3; i < 6; i++) { + if (set[i].style == FINAL) set[i].ftilt *= map[i]; + else if (set[i].style == DELTA) set[i].dtilt *= map[i]; + else if (set[i].style == VEL) set[i].vel *= map[i]; + } + + // set initial/final values for box size and shape for FINAL,DELTA,SCALE + // final not possible for VEL,RATE since don't know length of run yet + // final = initial if no setting + + for (int i = 0; i < 3; i++) { + set[i].lo_target = set[i].lo_stop = set[i].lo_start = domain->boxlo[i]; + set[i].hi_target = set[i].hi_stop = set[i].hi_start = domain->boxhi[i]; + + if (set[i].style == FINAL) { + set[i].lo_stop = set[i].flo; + set[i].hi_stop = set[i].fhi; + } else if (set[i].style == DELTA) { + set[i].lo_stop = set[i].lo_start + set[i].dlo; + set[i].hi_stop = set[i].hi_start + set[i].dhi; + } else if (set[i].style == SCALE) { + set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - + 0.5*set[i].scale*(set[i].hi_start-set[i].lo_start); + set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + + 0.5*set[i].scale*(set[i].hi_start-set[i].lo_start); + } + } + + for (int i = 3; i < 6; i++) { + if (i == 5) set[i].tilt_start = domain->xy; + else if (i == 4) set[i].tilt_start = domain->xz; + else if (i == 3) set[i].tilt_start = domain->yz; + set[i].tilt_flip = set[i].tilt_target = + set[i].tilt_stop = set[i].tilt_start; + + if (set[i].style == FINAL) { + set[i].tilt_stop = set[i].ftilt; + } else if (set[i].style == DELTA) { + set[i].tilt_stop = set[i].tilt_start + set[i].dtilt; + } + } + + // for VOLUME, setup links to other dims + // fixed, dynamic1,2, vol_start + + for (int i = 0; i < 3; i++) { + set[i].vol_start = domain->xprd * domain->yprd * domain->zprd; + + if (set[i].style != VOLUME) continue; + int other1 = (i+1) % 3; + int other2 = (i+2) % 3; + + if (set[other1].style == NONE) { + if (set[other2].style == NONE || set[other2].style == VOLUME) + error->all("Fix deform volume setting is invalid"); + set[i].substyle = ONE_FROM_ONE; + set[i].fixed = other1; + set[i].dynamic1 = other2; + } else if (set[other2].style == NONE) { + if (set[other1].style == NONE || set[other1].style == VOLUME) + error->all("Fix deform volume setting is invalid"); + set[i].substyle = ONE_FROM_ONE; + set[i].fixed = other2; + set[i].dynamic1 = other1; + } else if (set[other1].style == VOLUME) { + if (set[other2].style == NONE || set[other2].style == VOLUME) + error->all("Fix deform volume setting is invalid"); + set[i].substyle = TWO_FROM_ONE; + set[i].fixed = other1; + set[i].dynamic1 = other2; + } else if (set[other2].style == VOLUME) { + if (set[other1].style == NONE || set[other1].style == VOLUME) + error->all("Fix deform volume setting is invalid"); + set[i].substyle = TWO_FROM_ONE; + set[i].fixed = other2; + set[i].dynamic1 = other1; + } else { + set[i].substyle = ONE_FROM_TWO; + set[i].dynamic2 = other1; + set[i].dynamic2 = other2; + } + } + + // initial settings + // reneighboring only forced if flips will occur due to shape changes + + if (set[3].style || set[4].style || set[5].style) force_reneighbor = 1; + next_reneighbor = -1; + + nrigid = 0; + rfix = NULL; + flip = 0; +} + +/* ---------------------------------------------------------------------- */ + +FixDeform::~FixDeform() +{ + delete [] set; + delete [] rfix; +} + +/* ---------------------------------------------------------------------- */ + +int FixDeform::setmask() +{ + int mask = 0; + mask |= PRE_EXCHANGE; + mask |= END_OF_STEP; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixDeform::init() +{ + // error if more than one fix deform + // domain, fix nvt/sllod, compute temp/deform only work on single h_rate + + int count = 0; + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"deform") == 0) count++; + if (count > 1) error->warning("More than one fix deform"); + + // Kspace setting + + if (force->kspace) kspace_flag = 1; + else kspace_flag = 0; + + // elapsed time for entire simulation, including multiple runs if defined + // cannot be 0.0 since can't deform a finite distance in 0 time + + double delt = (update->endstep - update->beginstep) * update->dt; + if (delt == 0.0) error->all("Cannot use fix deform for 0 timestep run"); + + // set final values for box size and shape for VEL,RATE + // now possible since length of run is known + + for (int i = 0; i < 3; i++) { + if (set[i].style == VEL) { + set[i].lo_stop = set[i].lo_start - 0.5*delt*set[i].vel; + set[i].hi_stop = set[i].hi_start + 0.5*delt*set[i].vel; + } else if (set[i].style == RATE) { + set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - + 0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt)); + set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + + 0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt)); + } + } + for (int i = 3; i < 6; i++) { + if (set[i].style == VEL) { + set[i].tilt_stop = set[i].tilt_start + delt*set[i].vel; + } else if (set[i].style == RATE) { + set[i].tilt_stop = set[i].tilt_start * pow(1.0+set[i].rate,delt); + } + } + + // if using tilt RATE, then initial tilt must be non-zero + + for (int i = 3; i < 6; i++) + if (set[i].style == RATE) + if ((i == 5 && domain->xy == 0.0) || + (i == 4 && domain->xz == 0.0) || + (i == 3 && domain->yz == 0.0)) + error->all("Cannot use fix deform rate to tilt a box with zero tilt"); + + // if yz changes and will cause box flip, then xy cannot be changing + // this is b/c the flips would induce continuous changes in xz + // in order to keep the edge vectors of the flipped shape matrix + // a linear combination of the edge vectors of the unflipped shape matrix + + if (set[3].style && set[5].style) + if (set[3].tilt_stop < -0.5*(set[1].hi_start-set[1].lo_start) || + set[3].tilt_stop > 0.5*(set[1].hi_start-set[1].lo_start) || + set[3].tilt_stop < -0.5*(set[1].hi_stop-set[1].lo_stop) || + set[3].tilt_stop > 0.5*(set[1].hi_stop-set[1].lo_stop)) + error->all("Fix deform is changing yz by too much with changing xy"); + + // set domain->h_rate values for use by domain and other fixes/computes + // initialize all rates to 0.0 + // cannot set rate now for RATE,VOLUME styles since not constant + + h_rate = domain->h_rate; + h_ratelo = domain->h_ratelo; + + for (int i = 0; i < 3; i++) { + h_rate[i] = h_ratelo[i] = 0.0; + if (set[i].style == FINAL || set[i].style == DELTA || + set[i].style == SCALE || set[i].style == VEL) { + double dlo_dt = (set[i].lo_stop - set[i].lo_start) / delt; + double dhi_dt = (set[i].hi_stop - set[i].hi_start) / delt; + h_rate[i] = dhi_dt - dlo_dt; + h_ratelo[i] = dlo_dt; + } + } + for (int i = 3; i < 6; i++) { + h_rate[i] = 0.0; + if (set[i].style == FINAL || set[i].style == DELTA || + set[i].style == VEL) + h_rate[i] = (set[i].tilt_stop - set[i].tilt_start) / delt; + } + + // detect if any rigid fixes exist so rigid bodies can be rescaled + // rfix[] = indices to each fix rigid + + delete [] rfix; + nrigid = 0; + rfix = NULL; + + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) nrigid++; + if (nrigid) { + rfix = new int[nrigid]; + nrigid = 0; + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; + } +} + +/* ---------------------------------------------------------------------- + box flipped on previous step + perform irregular comm to migrate atoms to new procs + reset box tilts for flipped config and create new box in domain + remap to put far-away atoms back into new box + perform irregular on atoms in lamda coords to get atoms to new procs + force reneighboring on next timestep +------------------------------------------------------------------------- */ + +void FixDeform::pre_exchange() +{ + if (flip == 0) return; + + domain->yz = set[3].tilt_target = set[3].tilt_flip; + domain->xz = set[4].tilt_target = set[4].tilt_flip; + domain->xy = set[5].tilt_target = set[5].tilt_flip; + domain->set_global_box(); + domain->set_local_box(); + + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) + domain->remap(x[i],image[i]); + + domain->x2lamda(atom->nlocal); + comm->irregular(); + domain->lamda2x(atom->nlocal); + + flip = 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixDeform::end_of_step() +{ + int i; + + double delta = update->ntimestep - update->beginstep; + delta /= update->endstep - update->beginstep; + + // set new box size + // for RATE, set target directly based on current time and set h_rate + // for others except VOLUME, target is linear value between start and stop + + for (i = 0; i < 3; i++) { + if (set[i].style == NONE) continue; + + if (set[i].style == RATE) { + double delt = (update->ntimestep - update->beginstep) * update->dt; + set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - + 0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt)); + set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + + 0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt)); + h_rate[i] = set[i].rate * domain->h[i]; + h_ratelo[i] = -0.5*h_rate[i]; + + } else if (set[i].style != VOLUME) { + set[i].lo_target = set[i].lo_start + + delta*(set[i].lo_stop - set[i].lo_start); + set[i].hi_target = set[i].hi_start + + delta*(set[i].hi_stop - set[i].hi_start); + } + } + + // set new box size for VOLUME dims that are linked to other dims + // also need to set h_rate + + for (int i = 0; i < 3; i++) { + if (set[i].style != VOLUME) continue; + + if (set[i].substyle == ONE_FROM_ONE) { + set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - + 0.5*(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start)); + set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + + 0.5*(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start)); + + } else if (set[i].substyle == ONE_FROM_TWO) { + set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - + 0.5*(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].dynamic2].hi_target - + set[set[i].dynamic2].lo_target)); + set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + + 0.5*(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].dynamic2].hi_target - + set[set[i].dynamic2].lo_target)); + + } else if (set[i].substyle == TWO_FROM_ONE) { + set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - + 0.5*sqrt(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].fixed].hi_start - + set[set[i].fixed].lo_start) * + (set[i].hi_start - set[i].lo_start)); + set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + + 0.5*sqrt(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].fixed].hi_start - + set[set[i].fixed].lo_start) * + (set[i].hi_start - set[i].lo_start)); + } + } + + // for triclinic, set new box shape + // for RATE, set target directly based on current time and set h_rate + // for all others, target is linear value between start and stop values + + if (triclinic) { + double *h = domain->h; + + for (i = 3; i < 6; i++) { + if (set[i].style == NONE) continue; + + if (set[i].style == RATE) { + double delt = (update->ntimestep - update->beginstep) * update->dt; + set[i].tilt_target = set[i].tilt_start * pow(1.0+set[i].rate,delt); + h_rate[i] = set[i].rate * domain->h[i]; + } else if (set[i].style) { + set[i].tilt_target = set[i].tilt_start + + delta*(set[i].tilt_stop - set[i].tilt_start); + } + + // tilt_target can be large positive or large negative value + // add/subtract box lengths until tilt_target is closest to current value + + int idenom; + if (i == 5) idenom = 0; + else if (i == 4) idenom = 0; + else if (i == 3) idenom = 1; + double denom = set[idenom].hi_target - set[idenom].lo_target; + + double current = h[i]/h[idenom]; + while (set[i].tilt_target/denom - current > 0.0) + set[i].tilt_target -= denom; + while (set[i].tilt_target/denom - current < 0.0) + set[i].tilt_target += denom; + if (fabs(set[i].tilt_target/denom - 1.0 - current) < + fabs(set[i].tilt_target/denom - current)) + set[i].tilt_target -= denom; + } + } + + // if any tilt targets exceed bounds, set flip flag and new tilt_flip values + // flip will be performed on next timestep before reneighboring + // when yz flips and xy is non-zero, xz must also change + // this is to keep the edge vectors of the flipped shape matrix + // a linear combination of the edge vectors of the unflipped shape matrix + + if (triclinic) { + double xprd = set[0].hi_target - set[0].lo_target; + double yprd = set[1].hi_target - set[1].lo_target; + if (set[3].tilt_target < -0.5*yprd || set[3].tilt_target > 0.5*yprd || + set[4].tilt_target < -0.5*xprd || set[4].tilt_target > 0.5*xprd || + set[5].tilt_target < -0.5*xprd || set[5].tilt_target > 0.5*xprd) { + + flip = 1; + next_reneighbor = update->ntimestep + 1; + set[3].tilt_flip = set[3].tilt_target; + set[4].tilt_flip = set[4].tilt_target; + set[5].tilt_flip = set[5].tilt_target; + + if (set[3].tilt_flip < -0.5*yprd) { + set[3].tilt_flip += yprd; + set[4].tilt_flip += set[5].tilt_flip; + } else if (set[3].tilt_flip >= 0.5*yprd) { + set[3].tilt_flip -= yprd; + set[4].tilt_flip -= set[5].tilt_flip; + } + if (set[4].tilt_flip < -0.5*xprd) set[4].tilt_flip += xprd; + if (set[4].tilt_flip > 0.5*xprd) set[4].tilt_flip -= xprd; + if (set[5].tilt_flip < -0.5*xprd) set[5].tilt_flip += xprd; + if (set[5].tilt_flip > 0.5*xprd) set[5].tilt_flip -= xprd; + } + } + + // convert atoms to lamda coords + + if (remapflag == X_REMAP) { + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + domain->x2lamda(x[i],x[i]); + } + + // reset global and local box to new size/shape + + domain->boxlo[0] = set[0].lo_target; + domain->boxlo[1] = set[1].lo_target; + domain->boxlo[2] = set[2].lo_target; + domain->boxhi[0] = set[0].hi_target; + domain->boxhi[1] = set[1].hi_target; + domain->boxhi[2] = set[2].hi_target; + + if (triclinic) { + domain->yz = set[3].tilt_target; + domain->xz = set[4].tilt_target; + domain->xy = set[5].tilt_target; + } + + domain->set_global_box(); + domain->set_local_box(); + + // convert affine atoms back to box coords + + if (remapflag == X_REMAP) { + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + domain->lamda2x(x[i],x[i]); + } + + // remap centers-of-mass of rigid bodies + // needs to be new call, not dilate + // check where else fix dilate is called from + + /* + if (nrigid) + for (i = 0; i < nrigid; i++) + modify->fix[rfix[i]]->dilate(2,oldlo,oldhi,newlo,newhi); + */ + + // redo KSpace coeffs since box has changed + + if (kspace_flag) force->kspace->setup(); +} + +/* ---------------------------------------------------------------------- */ + +void FixDeform::options(int narg, char **arg) +{ + if (narg < 0) error->all("Illegal fix deform command"); + + remapflag = X_REMAP; + scaleflag = 1; + + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"remap") == 0) { + if (iarg+2 > narg) error->all("Illegal fix deform command"); + if (strcmp(arg[iarg+1],"x") == 0) remapflag = X_REMAP; + else if (strcmp(arg[iarg+1],"v") == 0) remapflag = V_REMAP; + else if (strcmp(arg[iarg+1],"none") == 0) remapflag = NO_REMAP; + else error->all("Illegal fix deform command"); + iarg += 2; + } else if (strcmp(arg[iarg],"units") == 0) { + if (iarg+2 > narg) error->all("Illegal fix deform command"); + if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; + else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; + else error->all("Illegal fix deform command"); + iarg += 2; + } else error->all("Illegal fix deform command"); + } +} diff --git a/src/fix_uniaxial.h b/src/fix_deform.h similarity index 53% rename from src/fix_uniaxial.h rename to src/fix_deform.h index c06f5add4b..bbaabbb4ea 100644 --- a/src/fix_uniaxial.h +++ b/src/fix_deform.h @@ -11,41 +11,54 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef FIX_UNIAXIAL_H -#define FIX_UNIAXIAL_H +#ifndef FIX_DEFORM_H +#define FIX_DEFORM_H #include "fix.h" namespace LAMMPS_NS { -class FixUniaxial : public Fix { +class FixDeform : public Fix { public: - FixUniaxial(class LAMMPS *, int, char **); - ~FixUniaxial(); + int remapflag; + + FixDeform(class LAMMPS *, int, char **); + ~FixDeform(); int setmask(); void init(); + void pre_exchange(); void end_of_step(); private: - int dir; - double lambda_final; - double alpha0; // asymmetry parameter - - double lambdai[3]; // initial - double lambdaf[3]; // final - double domainloi[3]; // initial box lo/hi - double domainhii[3]; - double *domainlo[3]; // pointers to make loop over dims possible - double *domainhi[3]; - double *domainprd[3]; - + double xlo_start,xhi_start,ylo_start,yhi_start,zlo_start,zhi_start; + double xlo_stop,xhi_stop,ylo_stop,yhi_stop,zlo_stop,zhi_stop; + double xy_start,xz_start,yz_start; + double xy_stop,xz_stop,yz_stop; + double xscale,yscale,zscale; + int triclinic,scaleflag,flip; + double *h_rate,*h_ratelo; + + struct Set { + int style,substyle; + double flo,fhi,ftilt; + double dlo,dhi,dtilt; + double scale,vel,rate; + double lo_start,hi_start; + double lo_stop,hi_stop; + double lo_target,hi_target; + double tilt_start,tilt_stop,tilt_target; + double vol_start,tilt_flip; + int fixed,dynamic1,dynamic2; + }; + Set *set; + int kspace_flag; // 1 if KSpace invoked, 0 if not int nrigid; // number of rigid fixes int *rfix; // indices of rigid fixes -}; -#endif + void options(int, char **); +}; } - +#endif diff --git a/src/fix_nvt_sllod.cpp b/src/fix_nvt_sllod.cpp new file mode 100644 index 0000000000..a18357a90d --- /dev/null +++ b/src/fix_nvt_sllod.cpp @@ -0,0 +1,278 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Pieter in 't Veld (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "stdlib.h" +#include "fix_nvt_sllod.h" +#include "math_extra.h" +#include "atom.h" +#include "force.h" +#include "domain.h" +#include "comm.h" +#include "group.h" +#include "update.h" +#include "respa.h" +#include "modify.h" +#include "fix.h" +#include "fix_deform.h" +#include "compute.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp + +/* ---------------------------------------------------------------------- */ + +FixNVTSlodd::FixNVTSlodd(LAMMPS *lmp, int narg, char **arg) : + FixNVT(lmp, narg, arg) {} + +/* ---------------------------------------------------------------------- */ + +void FixNVTSlodd::init() +{ + FixNVT::init(); + + // check fix deform remap settings + + int i; + for (i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"deform") == 0) { + if (((FixDeform *) modify->fix[i])->remapflag == X_REMAP && + comm->me == 0) + error->warning("Using compute temp/deform with fix deform remapping coords"); + break; + } + if (i == modify->nfix && comm->me == 0) + error->warning("Using compute temp/deform with no fix deform defined"); +} + +/* ---------------------------------------------------------------------- */ + +void FixNVTSlodd::initial_integrate() +{ + double dtfm; + + double delta = update->ntimestep - update->firststep; + delta /= update->nsteps; + t_target = t_start + delta * (t_stop-t_start); + + // update eta_dot + + f_eta = t_freq*t_freq * (t_current/t_target - 1.0); + eta_dot += f_eta*dthalf; + eta_dot *= drag_factor; + eta += dtv*eta_dot; + factor = exp(-dthalf*eta_dot); + + // update vthermal and x of only atoms in NVT group + // lamda = 0-1 triclinic lamda coords + // vstream = streaming velocity = Hrate*lamda + Hratelo + // vthermal = thermal velocity = v - vstream + // vdelu = Hrate*Hinv*vthermal + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *h_rate = domain->h_rate; + double *h_ratelo = domain->h_ratelo; + double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3]; + MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + vthermal[0] = v[i][0] - vstream[0]; + vthermal[1] = v[i][1] - vstream[1]; + vthermal[2] = v[i][2] - vstream[2]; + vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + + h_two[4]*vthermal[2]; + vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2]; + vdelu[2] = h_two[2]*vthermal[2]; + + v[i][0] = vstream[0] + + vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0]; + v[i][1] = vstream[1] + + vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1]; + v[i][2] = vstream[2] + + vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNVTSlodd::final_integrate() +{ + double dtfm; + + // update vthermal of only atoms in NVT group + // lamda = 0-1 triclinic lamda coords + // vstream = streaming velocity = Hrate*lamda + Hratelo + // vthermal = thermal velocity = v - vstream + // vdelu = Hrate*Hinv*vthermal + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *h_rate = domain->h_rate; + double *h_ratelo = domain->h_ratelo; + double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3]; + MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + vthermal[0] = v[i][0] - vstream[0]; + vthermal[1] = v[i][1] - vstream[1]; + vthermal[2] = v[i][2] - vstream[2]; + vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + + h_two[4]*vthermal[2]; + vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2]; + vdelu[2] = h_two[2]*vthermal[2]; + + v[i][0] = vstream[0] + + vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0]; + v[i][1] = vstream[1] + + vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1]; + v[i][2] = vstream[2] + + vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2]; + } + } + + // compute current T + + t_current = temperature->compute_scalar(); + + // update eta_dot + + f_eta = t_freq*t_freq * (t_current/t_target - 1.0); + eta_dot += f_eta*dthalf; + eta_dot *= drag_factor; +} + +/* ---------------------------------------------------------------------- */ + +void FixNVTSlodd::initial_integrate_respa(int ilevel, int flag) +{ + if (flag) return; // only used by NPT,NPH + + // set timesteps by level + + double dtfm; + dtv = step_respa[ilevel]; + dtf = 0.5 * step_respa[ilevel] * force->ftm2v; + dthalf = 0.5 * step_respa[ilevel]; + + // atom quantities + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + // outermost level - update eta_dot and apply to v + // all other levels - NVE update of v + + if (ilevel == nlevels_respa-1) { + double delta = update->ntimestep - update->firststep; + delta /= update->nsteps; + t_target = t_start + delta * (t_stop-t_start); + + // update eta_dot + + f_eta = t_freq*t_freq * (t_current/t_target - 1.0); + eta_dot += f_eta*dthalf; + eta_dot *= drag_factor; + eta += dtv*eta_dot; + factor = exp(-dthalf*eta_dot); + } else factor = 1.0; + + // update v of only atoms in NVT group + + double *h_rate = domain->h_rate; + double *h_ratelo = domain->h_ratelo; + double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3]; + MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + vthermal[0] = v[i][0] - vstream[0]; + vthermal[1] = v[i][1] - vstream[1]; + vthermal[2] = v[i][2] - vstream[2]; + vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + + h_two[4]*vthermal[2]; + vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2]; + vdelu[2] = h_two[2]*vthermal[2]; + + v[i][0] = vstream[0] + + vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0]; + v[i][1] = vstream[1] + + vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1]; + v[i][2] = vstream[2] + + vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2]; + } + } + + // innermost level - also update x of only atoms in NVT group + + if (ilevel == 0) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } + } +} diff --git a/src/fix_nvt_sllod.h b/src/fix_nvt_sllod.h new file mode 100644 index 0000000000..abf96f8380 --- /dev/null +++ b/src/fix_nvt_sllod.h @@ -0,0 +1,33 @@ +/* ---------------------------------------------------------------------- + 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_NVT_SLODD_H +#define FIX_NVT_SLODD_H + +#include "fix_nvt.h" + +namespace LAMMPS_NS { + +class FixNVTSlodd : public FixNVT { + public: + FixNVTSlodd(class LAMMPS *, int, char **); + ~FixNVTSlodd() {} + void init(); + void initial_integrate(); + void final_integrate(); + void initial_integrate_respa(int,int); +}; + +} + +#endif diff --git a/src/fix_uniaxial.cpp b/src/fix_uniaxial.cpp deleted file mode 100644 index 8b610fbeae..0000000000 --- a/src/fix_uniaxial.cpp +++ /dev/null @@ -1,212 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Carsten Svaneborg - (Max Planck Institute for Complex Systems, Dresden, Germany) -------------------------------------------------------------------------- */ - -#include "string.h" -#include "stdlib.h" -#include "math.h" -#include "fix_uniaxial.h" -#include "atom.h" -#include "update.h" -#include "force.h" -#include "domain.h" -#include "modify.h" -#include "comm.h" -#include "kspace.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -FixUniaxial::FixUniaxial(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg != 6) error->all("Illegal fix uniaxial command"); - - box_change = 1; - - nevery = atoi(arg[3]); - if (nevery <= 0) error->all("Illegal fix uniaxial command"); - - if (strcmp(arg[4],"x") == 0) dir = 0; - else if (strcmp(arg[4],"y") == 0) dir = 1; - else if (strcmp(arg[4],"z") == 0) dir = 2; - else error->all("Illegal fix uniaxial command"); - lambda_final = atof(arg[5]); - - if (lambda_final <= 0) error->all("Illegal fix uniaxial command"); - - if (domain->nonperiodic) - error->all("Cannot use fix uniaxial on non-periodic system"); - if (domain->triclinic) - error->all("Cannot use fix uniaxial with triclinic box"); - - nrigid = 0; - rfix = NULL; -} - -/* ---------------------------------------------------------------------- */ - -FixUniaxial::~FixUniaxial() -{ - delete [] rfix; -} - -/* ---------------------------------------------------------------------- */ - -int FixUniaxial::setmask() -{ - int mask = 0; - mask |= END_OF_STEP; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixUniaxial::init() -{ - // store pointers to domain variable so can loop over dimensions - - domainlo[0] = &domain->boxlo[0]; - domainlo[1] = &domain->boxlo[1]; - domainlo[2] = &domain->boxlo[2]; - - domainhi[0] = &domain->boxhi[0]; - domainhi[1] = &domain->boxhi[1]; - domainhi[2] = &domain->boxhi[2]; - - domainprd[0] = &domain->xprd; - domainprd[1] = &domain->yprd; - domainprd[2] = &domain->zprd; - - double L = pow((domain->boxhi[0]-domain->boxlo[0])* - (domain->boxhi[1]-domain->boxlo[1])* - (domain->boxhi[2]-domain->boxlo[2]) ,1.0/3.0); - - // save box sizes for coordinate rescaling - // calculate strains and asymmetry parameter - // alpha=lampdai[first]/lampbdai[second] for the two perp directions - - alpha0 = 1.0; - for (int m = 0; m < 3; m++) { - domainloi[m] = *domainlo[m]; - domainhii[m] = *domainhi[m]; - lambdai[m] = (*domainhi[m] - *domainlo[m])/L; - lambdaf[m] = ( m==dir ? lambda_final : 1.0/sqrt(lambda_final) ) ; - if (m != dir) { - alpha0*= lambdai[m]; - alpha0=1.0/alpha0; - } - } - - if (comm->me == 0) { - if (screen) { - fprintf(screen,"Initial strain = %g %g %g\n", - lambdai[0],lambdai[1],lambdai[2]); - fprintf(screen,"Target strain = %g %g %g\n", - lambdaf[0],lambdaf[1],lambdaf[2]); - } - if (logfile) { - fprintf(logfile,"Initial strain = %g %g %g\n", - lambdai[0],lambdai[1],lambdai[2]); - fprintf(logfile,"Target strain = %g %g %g\n", - lambdaf[0],lambdaf[1],lambdaf[2]); - } - } - - if (force->kspace) kspace_flag = 1; - else kspace_flag = 0; - - // detect if any fix rigid exist so rigid bodies can be re-scaled - // rfix[] = indices to each fix rigid - - delete [] rfix; - nrigid = 0; - rfix = NULL; - - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"rigid") == 0 || - strcmp(modify->fix[i]->style,"poems") == 0) nrigid++; - if (nrigid) { - rfix = new int[nrigid]; - nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"rigid") == 0 || - strcmp(modify->fix[i]->style,"poems") == 0) rfix[nrigid++] = i; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixUniaxial::end_of_step() -{ - int i,m; - double oldlo,oldhi,newlo,newhi,ratio; - - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - - double lvalue[3]; - - // linear interpolation of strain in specified direction - - lvalue[dir] = lambdai[dir]*(1.0-delta) + lambdaf[dir]*delta; - - // linear interpolation of asymmetry parameter in the perp direction - - double alpha = alpha0*(1-delta) + delta; - - // calculate strains perpendicular to dir - - for (m = 0; m < 3; m++) - if (m != dir) { - lvalue[m] = sqrt(alpha/lvalue[dir]); - alpha=1.0/alpha; - } - - // apply appropriate rescaling in each dimension - - double **x = atom->x; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - for (m = 0; m < 3; m++) { - oldlo = *domainlo[m]; - oldhi = *domainhi[m]; - - newlo = domainloi[m] * lvalue[m]/lambdai[m]; - newhi = domainhii[m] * lvalue[m]/lambdai[m]; - ratio = (newhi - newlo) / *domainprd[m]; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - x[i][m] = newlo + (x[i][m] - oldlo) * ratio; - - *domainlo[m] = newlo; - *domainhi[m] = newhi; - domain->prd[m] = *domainprd[m] = newhi - newlo; - - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->dilate(m,oldlo,oldhi,newlo,newhi); - } - - // redo KSpace coeffs since volume has changed - - if (kspace_flag) force->kspace->setup(); -} diff --git a/src/fix_volume_rescale.cpp b/src/fix_volume_rescale.cpp deleted file mode 100644 index 4861487bef..0000000000 --- a/src/fix_volume_rescale.cpp +++ /dev/null @@ -1,194 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "string.h" -#include "stdlib.h" -#include "math.h" -#include "fix_volume_rescale.h" -#include "update.h" -#include "force.h" -#include "modify.h" -#include "kspace.h" -#include "domain.h" -#include "atom.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -FixVolRescale::FixVolRescale(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg < 7) error->all("Illegal fix volume/rescale command"); - - box_change = 1; - - nevery = atoi(arg[3]); - if (nevery <= 0) error->all("Illegal fix volume/rescale command"); - - xflag = yflag = zflag = 0; - int iarg = 4; - while (iarg < narg) { - if (strcmp(arg[iarg],"x") == 0) { - if (iarg+3 > narg) error->all("Illegal fix volume/rescale command"); - if (domain->xperiodic == 0) - error->all("Cannot fix volume/rescale on a non-periodic boundary"); - xflag = 1; - xlo_stop = atof(arg[iarg+1]); - xhi_stop = atof(arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg],"y") == 0) { - if (iarg+3 > narg) error->all("Illegal fix volume/rescale command"); - if (domain->yperiodic == 0) - error->all("Cannot fix volume/rescale on a non-periodic boundary"); - yflag = 1; - ylo_stop = atof(arg[iarg+1]); - yhi_stop = atof(arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg],"z") == 0) { - if (iarg+3 > narg) error->all("Illegal fix volume/rescale command"); - if (domain->zperiodic == 0) - error->all("Cannot fix volume/rescale on a non-periodic boundary"); - zflag = 1; - zlo_stop = atof(arg[iarg+1]); - zhi_stop = atof(arg[iarg+2]); - iarg += 3; - } else error->all("Illegal fix volume/rescale command"); - } - - if (domain->triclinic) - error->all("Cannot use fix volume/rescale with triclinic box"); - - nrigid = 0; - rfix = NULL; -} - -/* ---------------------------------------------------------------------- */ - -FixVolRescale::~FixVolRescale() -{ - delete [] rfix; -} - -/* ---------------------------------------------------------------------- */ - -int FixVolRescale::setmask() -{ - int mask = 0; - mask |= END_OF_STEP; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixVolRescale::init() -{ - xlo_start = domain->boxlo[0]; - xhi_start = domain->boxhi[0]; - ylo_start = domain->boxlo[1]; - yhi_start = domain->boxhi[1]; - zlo_start = domain->boxlo[2]; - zhi_start = domain->boxhi[2]; - - if (force->kspace) kspace_flag = 1; - else kspace_flag = 0; - - // detect if any fix rigid exist so rigid bodies can be rescaled - // rfix[] = indices to each fix rigid - - delete [] rfix; - nrigid = 0; - rfix = NULL; - - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"rigid") == 0 || - strcmp(modify->fix[i]->style,"poems") == 0) nrigid++; - if (nrigid) { - rfix = new int[nrigid]; - nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"rigid") == 0 || - strcmp(modify->fix[i]->style,"poems") == 0) rfix[nrigid++] = i; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixVolRescale::end_of_step() -{ - int i; - double oldlo,oldhi,newlo,newhi,ratio; - - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - - double **x = atom->x; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - if (xflag) { - oldlo = domain->boxlo[0]; - oldhi = domain->boxhi[0]; - newlo = xlo_start + delta * (xlo_stop-xlo_start); - newhi = xhi_start + delta * (xhi_stop-xhi_start); - ratio = (newhi - newlo) / domain->xprd; - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - x[i][0] = newlo + (x[i][0] - oldlo) * ratio; - domain->boxlo[0] = newlo; - domain->boxhi[0] = newhi; - domain->prd[0] = domain->xprd = newhi - newlo; - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->dilate(0,oldlo,oldhi,newlo,newhi); - } - - if (yflag) { - oldlo = domain->boxlo[1]; - oldhi = domain->boxhi[1]; - newlo = ylo_start + delta * (ylo_stop-ylo_start); - newhi = yhi_start + delta * (yhi_stop-yhi_start); - ratio = (newhi - newlo) / domain->yprd; - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - x[i][1] = newlo + (x[i][1] - oldlo) * ratio; - domain->boxlo[1] = newlo; - domain->boxhi[1] = newhi; - domain->prd[1] = domain->yprd = newhi - newlo; - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->dilate(1,oldlo,oldhi,newlo,newhi); - } - - if (zflag) { - oldlo = domain->boxlo[2]; - oldhi = domain->boxhi[2]; - newlo = zlo_start + delta * (zlo_stop-zlo_start); - newhi = zhi_start + delta * (zhi_stop-zhi_start); - ratio = (newhi - newlo) / domain->zprd; - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - x[i][2] = newlo + (x[i][2] - oldlo) * ratio; - domain->boxlo[2] = newlo; - domain->boxhi[2] = newhi; - domain->prd[2] = domain->zprd = newhi - newlo; - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->dilate(2,oldlo,oldhi,newlo,newhi); - } - - // redo KSpace coeffs since volume has changed - - if (kspace_flag) force->kspace->setup(); -} diff --git a/src/math_extra.h b/src/math_extra.h new file mode 100755 index 0000000000..e043295681 --- /dev/null +++ b/src/math_extra.h @@ -0,0 +1,435 @@ +/* ---------------------------------------------------------------------- + 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 author: Mike Brown (SNL) +------------------------------------------------------------------------- */ + +#ifndef MATH_EXTRA_H +#define MATH_EXTRA_H + +#include "math.h" +#include "stdio.h" +#include "string.h" +#include "error.h" + +namespace MathExtra { + + // 3 vector operations + + inline void normalize3(const double *v, double *ans); + inline double dot3(const double *v1, const double *v2); + inline void cross3(const double *v1, const double *v2, double *ans); + + // 3x3 matrix operations + + inline double det3(const double mat[3][3]); + inline void diag_times3(const double *diagonal, const double mat[3][3], + double ans[3][3]); + inline void plus3(const double m[3][3], const double m2[3][3], + double ans[3][3]); + inline void transpose_times3(const double mat1[3][3], + const double mat2[3][3], + double ans[3][3]); + inline void invert3(const double mat[3][3], double ans[3][3]); + inline void times_column3(const double mat[3][3], const double*vec, + double *ans); + inline void transpose_times_column3(const double mat[3][3],const double*vec, + double *ans); + inline void row_times3(const double *v, const double m[3][3], + double *ans); + inline void mldivide3(const double mat[3][3], const double *vec, + double *ans, LAMMPS_NS::Error *error); + inline void write3(const double mat[3][3]); + + // quaternion operations + + inline void normalize4(double *quat); + inline void quat_to_mat(const double *quat, double mat[3][3]); + inline void quat_to_mat_trans(const double *quat, double mat[3][3]); + inline void axisangle_to_quat(const double *v, const double angle, + double *quat); + inline void multiply_quat_quat(const double *one, const double *two, + double *ans); + inline void multiply_vec_quat(const double *one, const double *two, + double *ans); + + // shape matrix operations + + inline void multiply_shape_shape(const double *one, const double *two, + double *ans); +}; + +/* ---------------------------------------------------------------------- + normalize a vector +------------------------------------------------------------------------- */ + +void MathExtra::normalize3(const double *v, double *ans) +{ + double den = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); + ans[0] = v[0]/den; + ans[1] = v[1]/den; + ans[2] = v[2]/den; +} + +/* ---------------------------------------------------------------------- + dot product of 2 vectors +------------------------------------------------------------------------- */ + +double MathExtra::dot3(const double *v1, const double *v2) +{ + return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]; +} + +/* ---------------------------------------------------------------------- + cross product of 2 vectors +------------------------------------------------------------------------- */ + +void MathExtra::cross3(const double *v1, const double *v2, double *ans) +{ + ans[0] = v1[1]*v2[2]-v1[2]*v2[1]; + ans[1] = v1[2]*v2[0]-v1[0]*v2[2]; + ans[2] = v1[0]*v2[1]-v1[1]*v2[0]; +} + +/* ---------------------------------------------------------------------- + determinant of a matrix +------------------------------------------------------------------------- */ + +double MathExtra::det3(const double m[3][3]) +{ + double ans = m[0][0]*m[1][1]*m[2][2] - m[0][0]*m[1][2]*m[2][1] - + m[1][0]*m[0][1]*m[2][2] + m[1][0]*m[0][2]*m[2][1] + + m[2][0]*m[0][1]*m[1][2] - m[2][0]*m[0][2]*m[1][1]; + return ans; +} + +/* ---------------------------------------------------------------------- + diagonal matrix times a full matrix +------------------------------------------------------------------------- */ + +void MathExtra::diag_times3(const double *d, const double m[3][3], + double ans[3][3]) +{ + ans[0][0] = d[0]*m[0][0]; + ans[0][1] = d[0]*m[0][1]; + ans[0][2] = d[0]*m[0][2]; + ans[1][0] = d[1]*m[1][0]; + ans[1][1] = d[1]*m[1][1]; + ans[1][2] = d[1]*m[1][2]; + ans[2][0] = d[2]*m[2][0]; + ans[2][1] = d[2]*m[2][1]; + ans[2][2] = d[2]*m[2][2]; +} + +/* ---------------------------------------------------------------------- + add two matrices +------------------------------------------------------------------------- */ + +void MathExtra::plus3(const double m[3][3], const double m2[3][3], + double ans[3][3]) +{ + ans[0][0] = m[0][0]+m2[0][0]; + ans[0][1] = m[0][1]+m2[0][1]; + ans[0][2] = m[0][2]+m2[0][2]; + ans[1][0] = m[1][0]+m2[1][0]; + ans[1][1] = m[1][1]+m2[1][1]; + ans[1][2] = m[1][2]+m2[1][2]; + ans[2][0] = m[2][0]+m2[2][0]; + ans[2][1] = m[2][1]+m2[2][1]; + ans[2][2] = m[2][2]+m2[2][2]; +} + +/* ---------------------------------------------------------------------- + multiply the transpose of mat1 times mat2 +------------------------------------------------------------------------- */ + +void MathExtra::transpose_times3(const double m[3][3], const double m2[3][3], + double ans[3][3]) +{ + ans[0][0] = m[0][0]*m2[0][0]+m[1][0]*m2[1][0]+m[2][0]*m2[2][0]; + ans[0][1] = m[0][0]*m2[0][1]+m[1][0]*m2[1][1]+m[2][0]*m2[2][1]; + ans[0][2] = m[0][0]*m2[0][2]+m[1][0]*m2[1][2]+m[2][0]*m2[2][2]; + ans[1][0] = m[0][1]*m2[0][0]+m[1][1]*m2[1][0]+m[2][1]*m2[2][0]; + ans[1][1] = m[0][1]*m2[0][1]+m[1][1]*m2[1][1]+m[2][1]*m2[2][1]; + ans[1][2] = m[0][1]*m2[0][2]+m[1][1]*m2[1][2]+m[2][1]*m2[2][2]; + ans[2][0] = m[0][2]*m2[0][0]+m[1][2]*m2[1][0]+m[2][2]*m2[2][0]; + ans[2][1] = m[0][2]*m2[0][1]+m[1][2]*m2[1][1]+m[2][2]*m2[2][1]; + ans[2][2] = m[0][2]*m2[0][2]+m[1][2]*m2[1][2]+m[2][2]*m2[2][2]; +} + +/* ---------------------------------------------------------------------- + invert a matrix + does NOT checks for singular or badly scaled matrix +------------------------------------------------------------------------- */ + +void MathExtra::invert3(const double m[3][3], double ans[3][3]) +{ + double den = m[0][0]*m[1][1]*m[2][2]-m[0][0]*m[1][2]*m[2][1]; + den += -m[1][0]*m[0][1]*m[2][2]+m[1][0]*m[0][2]*m[2][1]; + den += m[2][0]*m[0][1]*m[1][2]-m[2][0]*m[0][2]*m[1][1]; + + ans[0][0] = (m[1][1]*m[2][2]-m[1][2]*m[2][1]) / den; + ans[0][1] = -(m[0][1]*m[2][2]-m[0][2]*m[2][1]) / den; + ans[0][2] = (m[0][1]*m[1][2]-m[0][2]*m[1][1]) / den; + ans[1][0] = -(m[1][0]*m[2][2]-m[1][2]*m[2][0]) / den; + ans[1][1] = (m[0][0]*m[2][2]-m[0][2]*m[2][0]) / den; + ans[1][2] = -(m[0][0]*m[1][2]-m[0][2]*m[1][0]) / den; + ans[2][0] = (m[1][0]*m[2][1]-m[1][1]*m[2][0]) / den; + ans[2][1] = -(m[0][0]*m[2][1]-m[0][1]*m[2][0]) / den; + ans[2][2] = (m[0][0]*m[1][1]-m[0][1]*m[1][0]) / den; +} + +/* ---------------------------------------------------------------------- + matrix times vector +------------------------------------------------------------------------- */ + +void MathExtra::times_column3(const double m[3][3], const double *v, + double *ans) +{ + ans[0] = m[0][0]*v[0]+m[0][1]*v[1]+m[0][2]*v[2]; + ans[1] = m[1][0]*v[0]+m[1][1]*v[1]+m[1][2]*v[2]; + ans[2] = m[2][0]*v[0]+m[2][1]*v[1]+m[2][2]*v[2]; +} + +/* ---------------------------------------------------------------------- + transposed matrix times vector +------------------------------------------------------------------------- */ + +void MathExtra::transpose_times_column3(const double m[3][3], const double *v, + double *ans) +{ + ans[0] = m[0][0]*v[0]+m[1][0]*v[1]+m[2][0]*v[2]; + ans[1] = m[0][1]*v[0]+m[1][1]*v[1]+m[2][1]*v[2]; + ans[2] = m[0][2]*v[0]+m[1][2]*v[1]+m[2][2]*v[2]; +} + +/* ---------------------------------------------------------------------- + row vector times matrix +------------------------------------------------------------------------- */ + +void MathExtra::row_times3(const double *v, const double m[3][3], + double *ans) +{ + ans[0] = m[0][0]*v[0]+v[1]*m[1][0]+v[2]*m[2][0]; + ans[1] = v[0]*m[0][1]+m[1][1]*v[1]+v[2]*m[2][1]; + ans[2] = v[0]*m[0][2]+v[1]*m[1][2]+m[2][2]*v[2]; +} + +/* ---------------------------------------------------------------------- + solve Ax = b or M ans = v + use gaussian elimination & partial pivoting on matrix +------------------------------------------------------------------------- */ + +void MathExtra::mldivide3(const double m[3][3], const double *v, double *ans, + LAMMPS_NS::Error *error) +{ + // create augmented matrix for pivoting + + double aug[3][4]; + for (unsigned i = 0; i < 3; i++) { + aug[i][3] = v[i]; + for (unsigned j = 0; j < 3; j++) aug[i][j] = m[i][j]; + } + + for (unsigned i = 0; i < 2; i++) { + unsigned p = i; + for (unsigned j = i+1; j < 3; j++) { + if (fabs(aug[j][i]) > fabs(aug[i][i])) { + double tempv[4]; + memcpy(tempv,aug[i],4*sizeof(double)); + memcpy(aug[i],aug[j],4*sizeof(double)); + memcpy(aug[j],tempv,4*sizeof(double)); + } + } + + while (aug[p][i] == 0 && p < 3) p++; + + if (p == 3) error->all("Bad matrix inversion in MathExtra::mldivide3"); + else + if (p != i) { + double tempv[4]; + memcpy(tempv,aug[i],4*sizeof(double)); + memcpy(aug[i],aug[p],4*sizeof(double)); + memcpy(aug[p],tempv,4*sizeof(double)); + } + + for (unsigned j = i+1; j < 3; j++) { + double m = aug[j][i]/aug[i][i]; + for (unsigned k=i+1; k<4; k++) aug[j][k]-=m*aug[i][k]; + } + } + + if (aug[2][2] == 0) + error->all("Bad matrix inversion in MathExtra::mldivide3"); + + // back substitution + + ans[2] = aug[2][3]/aug[2][2]; + for (int i = 1; i >= 0; i--) { + double sumax = 0.0; + for (unsigned j = i+1; j < 3; j++) sumax += aug[i][j]*ans[j]; + ans[i] = (aug[i][3]-sumax) / aug[i][i]; + } +} + +/* ---------------------------------------------------------------------- + output a matrix +------------------------------------------------------------------------- */ + +void MathExtra::write3(const double mat[3][3]) +{ + for (unsigned i = 0; i < 3; i++) { + for (unsigned j = 0; j < 3; j++) printf("%g ",mat[i][j]); + printf("\n"); + } +} + +/* ---------------------------------------------------------------------- + normalize a quaternion +------------------------------------------------------------------------- */ + +void MathExtra::normalize4(double *quat) +{ + double den = sqrt(quat[0]*quat[0]+quat[1]*quat[1] + + quat[2]*quat[2]+quat[3]*quat[3]); + quat[0] /= den; + quat[1] /= den; + quat[2] /= den; + quat[3] /= den; +} + +/* ---------------------------------------------------------------------- + compute rotation matrix from quaternion + quat = [w i j k] +------------------------------------------------------------------------- */ + +void MathExtra::quat_to_mat(const double *quat, double mat[3][3]) +{ + double w2 = quat[0]*quat[0]; + double i2 = quat[1]*quat[1]; + double j2 = quat[2]*quat[2]; + double k2 = quat[3]*quat[3]; + double twoij = 2.0*quat[1]*quat[2]; + double twoik = 2.0*quat[1]*quat[3]; + double twojk = 2.0*quat[2]*quat[3]; + double twoiw = 2.0*quat[1]*quat[0]; + double twojw = 2.0*quat[2]*quat[0]; + double twokw = 2.0*quat[3]*quat[0]; + + mat[0][0] = w2+i2-j2-k2; + mat[0][1] = twoij-twokw; + mat[0][2] = twojw+twoik; + + mat[1][0] = twoij+twokw; + mat[1][1] = w2-i2+j2-k2; + mat[1][2] = twojk-twoiw; + + mat[2][0] = twoik-twojw; + mat[2][1] = twojk+twoiw; + mat[2][2] = w2-i2-j2+k2; +} + +/* ---------------------------------------------------------------------- + compute rotation matrix from quaternion conjugate + quat = [w i j k] +------------------------------------------------------------------------- */ + +void MathExtra::quat_to_mat_trans(const double *quat, double mat[3][3]) +{ + double w2 = quat[0]*quat[0]; + double i2 = quat[1]*quat[1]; + double j2 = quat[2]*quat[2]; + double k2 = quat[3]*quat[3]; + double twoij = 2.0*quat[1]*quat[2]; + double twoik = 2.0*quat[1]*quat[3]; + double twojk = 2.0*quat[2]*quat[3]; + double twoiw = 2.0*quat[1]*quat[0]; + double twojw = 2.0*quat[2]*quat[0]; + double twokw = 2.0*quat[3]*quat[0]; + + mat[0][0] = w2+i2-j2-k2; + mat[1][0] = twoij-twokw; + mat[2][0] = twojw+twoik; + + mat[0][1] = twoij+twokw; + mat[1][1] = w2-i2+j2-k2; + mat[2][1] = twojk-twoiw; + + mat[0][2] = twoik-twojw; + mat[1][2] = twojk+twoiw; + mat[2][2] = w2-i2-j2+k2; +} + +/* ---------------------------------------------------------------------- + compute quaternion from axis-angle rotation + v MUST be a unit vector +------------------------------------------------------------------------- */ + +void MathExtra::axisangle_to_quat(const double *v, const double angle, + double *quat) +{ + double halfa = 0.5*angle; + double sina = sin(halfa); + quat[0] = cos(halfa); + quat[1] = v[0]*sina; + quat[2] = v[1]*sina; + quat[3] = v[2]*sina; +} + +/* ---------------------------------------------------------------------- + multiply 2 quaternions + effectively concatenates rotations + NOT a commutative operation +------------------------------------------------------------------------- */ + +void MathExtra::multiply_quat_quat(const double *one, const double *two, + double *ans) +{ + ans[0] = one[0]*two[0]-one[1]*two[1]-one[2]*two[2]-one[3]*two[3]; + ans[1] = one[0]*two[1]+one[1]*two[0]+one[2]*two[3]-one[3]*two[2]; + ans[2] = one[0]*two[2]-one[1]*two[3]+one[2]*two[0]+one[3]*two[1]; + ans[3] = one[0]*two[3]+one[1]*two[2]-one[2]*two[1]+one[3]*two[0]; +} + +/* ---------------------------------------------------------------------- + multiply 3 vector times quaternion + 3 vector one is treated as quaternion (0,one) +------------------------------------------------------------------------- */ + +void MathExtra::multiply_vec_quat(const double *one, const double *two, + double *ans) +{ + ans[0] = -one[0]*two[1]-one[1]*two[2]-one[2]*two[3]; + ans[1] = one[0]*two[0]+one[1]*two[3]-one[2]*two[2]; + ans[2] = -one[0]*two[3]+one[1]*two[0]+one[2]*two[1]; + ans[3] = one[0]*two[2]-one[1]*two[1]+one[2]*two[0]; +} + +/* ---------------------------------------------------------------------- + multiply 2 shape matrices + upper-triangular 3x3, stored as 6-vector in Voigt notation +------------------------------------------------------------------------- */ + +void MathExtra::multiply_shape_shape(const double *one, const double *two, + double *ans) +{ + ans[0] = one[0]*two[0]; + ans[1] = one[1]*two[1]; + ans[2] = one[2]*two[2]; + ans[3] = one[1]*two[3] + one[3]*two[2]; + ans[4] = one[0]*two[4] + one[5]*two[3] + one[4]*two[2]; + ans[5] = one[0]*two[5] + one[5]*two[1]; +} + +#endif diff --git a/src/style_asphere.h b/src/style_asphere.h new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/style_class2.h b/src/style_class2.h index 3ba9b32df4..e69de29bb2 100644 --- a/src/style_class2.h +++ b/src/style_class2.h @@ -1,56 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef AngleInclude -#include "angle_class2.h" -#endif - -#ifdef AngleClass -AngleStyle(class2,AngleClass2) -#endif - -#ifdef BondInclude -#include "bond_class2.h" -#endif - -#ifdef BondClass -BondStyle(class2,BondClass2) -#endif - -#ifdef DihedralInclude -#include "dihedral_class2.h" -#endif - -#ifdef DihedralClass -DihedralStyle(class2,DihedralClass2) -#endif - -#ifdef ImproperInclude -#include "improper_class2.h" -#endif - -#ifdef ImproperClass -ImproperStyle(class2,ImproperClass2) -#endif - -#ifdef PairInclude -#include "pair_lj_class2.h" -#include "pair_lj_class2_coul_cut.h" -#include "pair_lj_class2_coul_long.h" -#endif - -#ifdef PairClass -PairStyle(lj/class2,PairLJClass2) -PairStyle(lj/class2/coul/cut,PairLJClass2CoulCut) -PairStyle(lj/class2/coul/long,PairLJClass2CoulLong) -#endif diff --git a/src/style_colloid.h b/src/style_colloid.h new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/style_dipole.h b/src/style_dipole.h new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/style_dpd.h b/src/style_dpd.h index 8ce617c0c2..e69de29bb2 100644 --- a/src/style_dpd.h +++ b/src/style_dpd.h @@ -1,28 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef AtomInclude -#include "atom_vec_dpd.h" -#endif - -#ifdef AtomClass -AtomStyle(dpd,AtomVecDPD) -#endif - -#ifdef PairInclude -#include "pair_dpd.h" -#endif - -#ifdef PairClass -PairStyle(dpd,PairDPD) -#endif diff --git a/src/style_granular.h b/src/style_granular.h index 7b0f4b6a71..e69de29bb2 100644 --- a/src/style_granular.h +++ b/src/style_granular.h @@ -1,50 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef AtomInclude -#include "atom_vec_granular.h" -#endif - -#ifdef AtomClass -AtomStyle(granular,AtomVecGranular) -# endif - -#ifdef FixInclude -#include "fix_freeze.h" -#include "fix_gran_diag.h" -#include "fix_nve_gran.h" -#include "fix_pour.h" -#include "fix_shear_history.h" -#include "fix_wall_gran.h" -#endif - -#ifdef FixClass -FixStyle(freeze,FixFreeze) -FixStyle(gran/diag,FixGranDiag) -FixStyle(nve/gran,FixNVEGran) -FixStyle(pour,FixPour) -FixStyle(SHEAR_HISTORY,FixShearHistory) -FixStyle(wall/gran,FixWallGran) -#endif - -#ifdef PairInclude -#include "pair_gran_hertzian.h" -#include "pair_gran_history.h" -#include "pair_gran_no_history.h" -#endif - -#ifdef PairClass -PairStyle(gran/hertzian,PairGranHertzian) -PairStyle(gran/history,PairGranHistory) -PairStyle(gran/no_history,PairGranNoHistory) -#endif diff --git a/src/style_kspace.h b/src/style_kspace.h index 52cae87bdb..e69de29bb2 100644 --- a/src/style_kspace.h +++ b/src/style_kspace.h @@ -1,38 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef KSpaceInclude -#include "ewald.h" -#include "pppm.h" -#include "pppm_tip4p.h" -#endif - -#ifdef KSpaceClass -KSpaceStyle(ewald,Ewald) -KSpaceStyle(pppm,PPPM) -KSpaceStyle(pppm/tip4p,PPPMTIP4P) -#endif - -#ifdef PairInclude -#include "pair_buck_coul_long.h" -#include "pair_lj_cut_coul_long.h" -#include "pair_lj_cut_coul_long_tip4p.h" -#include "pair_lj_charmm_coul_long.h" -#endif - -#ifdef PairClass -PairStyle(buck/coul/long,PairBuckCoulLong) -PairStyle(lj/cut/coul/long,PairLJCutCoulLong) -PairStyle(lj/cut/coul/long/tip4p,PairLJCutCoulLongTIP4P) -PairStyle(lj/charmm/coul/long,PairLJCharmmCoulLong) -#endif diff --git a/src/style_manybody.h b/src/style_manybody.h index c564993e3e..e69de29bb2 100644 --- a/src/style_manybody.h +++ b/src/style_manybody.h @@ -1,30 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef PairInclude -#include "pair_airebo.h" -#include "pair_eam.h" -#include "pair_eam_alloy.h" -#include "pair_eam_fs.h" -#include "pair_sw.h" -#include "pair_tersoff.h" -#endif - -#ifdef PairClass -PairStyle(airebo,PairAIREBO) -PairStyle(eam,PairEAM) -PairStyle(eam/alloy,PairEAMAlloy) -PairStyle(eam/fs,PairEAMFS) -PairStyle(sw,PairSW) -PairStyle(tersoff,PairTersoff) -#endif diff --git a/src/style_molecule.h b/src/style_molecule.h index 6e424035a9..e69de29bb2 100644 --- a/src/style_molecule.h +++ b/src/style_molecule.h @@ -1,116 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef AngleInclude -#include "angle_charmm.h" -#include "angle_cosine.h" -#include "angle_cosine_squared.h" -#include "angle_harmonic.h" -#include "angle_hybrid.h" -#endif - -#ifdef AngleClass -AngleStyle(charmm,AngleCharmm) -AngleStyle(cosine,AngleCosine) -AngleStyle(cosine/squared,AngleCosineSquared) -AngleStyle(harmonic,AngleHarmonic) -AngleStyle(hybrid,AngleHybrid) -#endif - -#ifdef AtomInclude -#include "atom_vec_angle.h" -#include "atom_vec_bond.h" -#include "atom_vec_full.h" -#include "atom_vec_molecular.h" -#endif - -#ifdef AtomClass -AtomStyle(angle,AtomVecAngle) -AtomStyle(bond,AtomVecBond) -AtomStyle(full,AtomVecFull) -AtomStyle(molecular,AtomVecMolecular) -#endif - -#ifdef BondInclude -#include "bond_fene.h" -#include "bond_fene_expand.h" -#include "bond_harmonic.h" -#include "bond_hybrid.h" -#include "bond_morse.h" -#include "bond_nonlinear.h" -#include "bond_quartic.h" -#endif - -#ifdef BondClass -BondStyle(fene,BondFENE) -BondStyle(fene/expand,BondFENEExpand) -BondStyle(harmonic,BondHarmonic) -BondStyle(hybrid,BondHybrid) -BondStyle(morse,BondMorse) -BondStyle(nonlinear,BondNonlinear) -BondStyle(quartic,BondQuartic) -#endif - -#ifdef DihedralInclude -#include "dihedral_charmm.h" -#include "dihedral_harmonic.h" -#include "dihedral_helix.h" -#include "dihedral_hybrid.h" -#include "dihedral_multi_harmonic.h" -#include "dihedral_opls.h" -#endif - -#ifdef DihedralClass -DihedralStyle(charmm,DihedralCharmm) -DihedralStyle(harmonic,DihedralHarmonic) -DihedralStyle(helix,DihedralHelix) -DihedralStyle(hybrid,DihedralHybrid) -DihedralStyle(multi/harmonic,DihedralMultiHarmonic) -DihedralStyle(opls,DihedralOPLS) -#endif - -#ifdef DumpInclude -#include "dump_bond.h" -#endif - -#ifdef DumpClass -DumpStyle(bond,DumpBond) -#endif - -#ifdef FixInclude -#endif - -#ifdef FixClass -#endif - -#ifdef ImproperInclude -#include "improper_cvff.h" -#include "improper_harmonic.h" -#include "improper_hybrid.h" -#endif - -#ifdef ImproperClass -ImproperStyle(cvff,ImproperCvff) -ImproperStyle(harmonic,ImproperHarmonic) -ImproperStyle(hybrid,ImproperHybrid) -#endif - -#ifdef PairInclude -#include "pair_lj_charmm_coul_charmm.h" -#include "pair_lj_charmm_coul_charmm_implicit.h" -#endif - -#ifdef PairClass -PairStyle(lj/charmm/coul/charmm,PairLJCharmmCoulCharmm) -PairStyle(lj/charmm/coul/charmm/implicit,PairLJCharmmCoulCharmmImplicit) -#endif diff --git a/src/style_opt.h b/src/style_opt.h index 061816e136..e69de29bb2 100644 --- a/src/style_opt.h +++ b/src/style_opt.h @@ -1,30 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef PairInclude -#include "pair_eam_opt.h" -#include "pair_eam_alloy_opt.h" -#include "pair_eam_fs_opt.h" -#include "pair_lj_charmm_coul_long_opt.h" -#include "pair_lj_cut_opt.h" -#include "pair_morse_opt.h" -#endif - -#ifdef PairClass -PairStyle(eam/opt,PairEAMOpt) -PairStyle(eam/alloy/opt,PairEAMAlloyOpt) -PairStyle(eam/fs/opt,PairEAMFSOpt) -PairStyle(lj/cut/opt,PairLJCutOpt) -PairStyle(lj/charmm/coul/long/opt,PairLJCharmmCoulLongOpt) -PairStyle(morse/opt,PairMorseOpt) -#endif diff --git a/src/style_xtc.h b/src/style_xtc.h index 7110dda312..e69de29bb2 100644 --- a/src/style_xtc.h +++ b/src/style_xtc.h @@ -1,20 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef DumpInclude -#include "dump_xtc.h" -#endif - -#ifdef DumpClass -DumpStyle(xtc,DumpXTC) -#endif