From 4cadf1cf4969987b70b50279b181d136969fe661 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 19 Jan 2010 00:04:52 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3727 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/angle.cpp | 222 ++++++++++++++++++++++++++++++++++++++++++ src/bond.cpp | 204 ++++++++++++++++++++++++++++++++++++++ src/dihedral.cpp | 249 +++++++++++++++++++++++++++++++++++++++++++++++ src/improper.cpp | 243 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 918 insertions(+) create mode 100644 src/angle.cpp create mode 100644 src/bond.cpp create mode 100644 src/dihedral.cpp create mode 100644 src/improper.cpp diff --git a/src/angle.cpp b/src/angle.cpp new file mode 100644 index 0000000000..df8a80e22f --- /dev/null +++ b/src/angle.cpp @@ -0,0 +1,222 @@ +/* ---------------------------------------------------------------------- + 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 "math.h" +#include "angle.h" +#include "atom.h" +#include "force.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +Angle::Angle(LAMMPS *lmp) : Pointers(lmp) +{ + energy = 0.0; + + allocated = 0; + PI = 4.0*atan(1.0); + THIRD = 1.0/3.0; + + maxeatom = maxvatom = 0; + eatom = NULL; + vatom = NULL; +} + +/* ---------------------------------------------------------------------- */ + +Angle::~Angle() +{ + memory->sfree(eatom); + memory->destroy_2d_double_array(vatom); +} + +/* ---------------------------------------------------------------------- + check if all coeffs are set +------------------------------------------------------------------------- */ + +void Angle::init() +{ + if (!allocated) error->all("Angle coeffs are not set"); + for (int i = 1; i <= atom->nangletypes; i++) + if (setflag[i] == 0) error->all("All angle coeffs are not set"); +} + +/* ---------------------------------------------------------------------- + setup for energy, virial computation + see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) +------------------------------------------------------------------------- */ + +void Angle::ev_setup(int eflag, int vflag) +{ + int i,n; + + evflag = 1; + + eflag_either = eflag; + eflag_global = eflag % 2; + eflag_atom = eflag / 2; + + vflag_either = vflag; + vflag_global = vflag % 4; + vflag_atom = vflag / 4; + + // reallocate per-atom arrays if necessary + + if (eflag_atom && atom->nmax > maxeatom) { + maxeatom = atom->nmax; + memory->sfree(eatom); + eatom = (double *) memory->smalloc(maxeatom*sizeof(double),"bond:eatom"); + } + if (vflag_atom && atom->nmax > maxvatom) { + maxvatom = atom->nmax; + memory->destroy_2d_double_array(vatom); + vatom = memory->create_2d_double_array(maxvatom,6,"bond:vatom"); + } + + // zero accumulators + + if (eflag_global) energy = 0.0; + if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; + if (eflag_atom) { + n = atom->nlocal; + if (force->newton_bond) n += atom->nghost; + for (i = 0; i < n; i++) eatom[i] = 0.0; + } + if (vflag_atom) { + n = atom->nlocal; + if (force->newton_bond) n += atom->nghost; + for (i = 0; i < n; i++) { + vatom[i][0] = 0.0; + vatom[i][1] = 0.0; + vatom[i][2] = 0.0; + vatom[i][3] = 0.0; + vatom[i][4] = 0.0; + vatom[i][5] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + tally energy and virial into global and per-atom accumulators + virial = r1F1 + r2F2 + r3F3 = (r1-r2) F1 + (r3-r2) F3 = del1*f1 + del2*f3 +------------------------------------------------------------------------- */ + +void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond, + double eangle, double *f1, double *f3, + double delx1, double dely1, double delz1, + double delx2, double dely2, double delz2) +{ + double eanglethird,v[6]; + + if (eflag_either) { + if (eflag_global) { + if (newton_bond) energy += eangle; + else { + eanglethird = THIRD*eangle; + if (i < nlocal) energy += eanglethird; + if (j < nlocal) energy += eanglethird; + if (k < nlocal) energy += eanglethird; + } + } + if (eflag_atom) { + eanglethird = THIRD*eangle; + if (newton_bond || i < nlocal) eatom[i] += eanglethird; + if (newton_bond || j < nlocal) eatom[j] += eanglethird; + if (newton_bond || k < nlocal) eatom[k] += eanglethird; + } + } + + if (vflag_either) { + v[0] = delx1*f1[0] + delx2*f3[0]; + v[1] = dely1*f1[1] + dely2*f3[1]; + v[2] = delz1*f1[2] + delz2*f3[2]; + v[3] = delx1*f1[1] + delx2*f3[1]; + v[4] = delx1*f1[2] + delx2*f3[2]; + v[5] = dely1*f1[2] + dely2*f3[2]; + + if (vflag_global) { + if (newton_bond) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } else { + if (i < nlocal) { + virial[0] += THIRD*v[0]; + virial[1] += THIRD*v[1]; + virial[2] += THIRD*v[2]; + virial[3] += THIRD*v[3]; + virial[4] += THIRD*v[4]; + virial[5] += THIRD*v[5]; + } + if (j < nlocal) { + virial[0] += THIRD*v[0]; + virial[1] += THIRD*v[1]; + virial[2] += THIRD*v[2]; + virial[3] += THIRD*v[3]; + virial[4] += THIRD*v[4]; + virial[5] += THIRD*v[5]; + } + if (k < nlocal) { + virial[0] += THIRD*v[0]; + virial[1] += THIRD*v[1]; + virial[2] += THIRD*v[2]; + virial[3] += THIRD*v[3]; + virial[4] += THIRD*v[4]; + virial[5] += THIRD*v[5]; + } + } + } + + if (vflag_atom) { + if (newton_bond || i < nlocal) { + vatom[i][0] += THIRD*v[0]; + vatom[i][1] += THIRD*v[1]; + vatom[i][2] += THIRD*v[2]; + vatom[i][3] += THIRD*v[3]; + vatom[i][4] += THIRD*v[4]; + vatom[i][5] += THIRD*v[5]; + } + if (newton_bond || j < nlocal) { + vatom[j][0] += THIRD*v[0]; + vatom[j][1] += THIRD*v[1]; + vatom[j][2] += THIRD*v[2]; + vatom[j][3] += THIRD*v[3]; + vatom[j][4] += THIRD*v[4]; + vatom[j][5] += THIRD*v[5]; + } + if (newton_bond || k < nlocal) { + vatom[k][0] += THIRD*v[0]; + vatom[k][1] += THIRD*v[1]; + vatom[k][2] += THIRD*v[2]; + vatom[k][3] += THIRD*v[3]; + vatom[k][4] += THIRD*v[4]; + vatom[k][5] += THIRD*v[5]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +double Angle::memory_usage() +{ + double bytes = maxeatom * sizeof(double); + bytes += maxvatom*6 * sizeof(double); + return bytes; +} diff --git a/src/bond.cpp b/src/bond.cpp new file mode 100644 index 0000000000..6551297f46 --- /dev/null +++ b/src/bond.cpp @@ -0,0 +1,204 @@ +/* ---------------------------------------------------------------------- + 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 "bond.h" +#include "atom.h" +#include "force.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ----------------------------------------------------------------------- + set bond contribution to Vdwl energy to 0.0 + a particular bond style can override this +------------------------------------------------------------------------- */ + +Bond::Bond(LAMMPS *lmp) : Pointers(lmp) +{ + energy = 0.0; + + allocated = 0; + + maxeatom = maxvatom = 0; + eatom = NULL; + vatom = NULL; +} + +/* ---------------------------------------------------------------------- */ + +Bond::~Bond() +{ + memory->sfree(eatom); + memory->destroy_2d_double_array(vatom); +} + +/* ---------------------------------------------------------------------- + check if all coeffs are set +------------------------------------------------------------------------- */ + +void Bond::init() +{ + if (!allocated) error->all("Bond coeffs are not set"); + for (int i = 1; i <= atom->nbondtypes; i++) + if (setflag[i] == 0) error->all("All bond coeffs are not set"); + init_style(); +} + +/* ---------------------------------------------------------------------- + setup for energy, virial computation + see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) +------------------------------------------------------------------------- */ + +void Bond::ev_setup(int eflag, int vflag) +{ + int i,n; + + evflag = 1; + + eflag_either = eflag; + eflag_global = eflag % 2; + eflag_atom = eflag / 2; + + vflag_either = vflag; + vflag_global = vflag % 4; + vflag_atom = vflag / 4; + + // reallocate per-atom arrays if necessary + + if (eflag_atom && atom->nmax > maxeatom) { + maxeatom = atom->nmax; + memory->sfree(eatom); + eatom = (double *) memory->smalloc(maxeatom*sizeof(double),"bond:eatom"); + } + if (vflag_atom && atom->nmax > maxvatom) { + maxvatom = atom->nmax; + memory->destroy_2d_double_array(vatom); + vatom = memory->create_2d_double_array(maxvatom,6,"bond:vatom"); + } + + // zero accumulators + + if (eflag_global) energy = 0.0; + if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; + if (eflag_atom) { + n = atom->nlocal; + if (force->newton_bond) n += atom->nghost; + for (i = 0; i < n; i++) eatom[i] = 0.0; + } + if (vflag_atom) { + n = atom->nlocal; + if (force->newton_bond) n += atom->nghost; + for (i = 0; i < n; i++) { + vatom[i][0] = 0.0; + vatom[i][1] = 0.0; + vatom[i][2] = 0.0; + vatom[i][3] = 0.0; + vatom[i][4] = 0.0; + vatom[i][5] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + tally energy and virial into global and per-atom accumulators +------------------------------------------------------------------------- */ + +void Bond::ev_tally(int i, int j, int nlocal, int newton_bond, + double ebond, double fbond, + double delx, double dely, double delz) +{ + double ebondhalf,v[6]; + + if (eflag_either) { + if (eflag_global) { + if (newton_bond) energy += ebond; + else { + ebondhalf = 0.5*ebond; + if (i < nlocal) energy += ebondhalf; + if (j < nlocal) energy += ebondhalf; + } + } + if (eflag_atom) { + ebondhalf = 0.5*ebond; + if (newton_bond || i < nlocal) eatom[i] += ebondhalf; + if (newton_bond || j < nlocal) eatom[j] += ebondhalf; + } + } + + if (vflag_either) { + v[0] = delx*delx*fbond; + v[1] = dely*dely*fbond; + v[2] = delz*delz*fbond; + v[3] = delx*dely*fbond; + v[4] = delx*delz*fbond; + v[5] = dely*delz*fbond; + + if (vflag_global) { + if (newton_bond) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } else { + if (i < nlocal) { + virial[0] += 0.5*v[0]; + virial[1] += 0.5*v[1]; + virial[2] += 0.5*v[2]; + virial[3] += 0.5*v[3]; + virial[4] += 0.5*v[4]; + virial[5] += 0.5*v[5]; + } + if (j < nlocal) { + virial[0] += 0.5*v[0]; + virial[1] += 0.5*v[1]; + virial[2] += 0.5*v[2]; + virial[3] += 0.5*v[3]; + virial[4] += 0.5*v[4]; + virial[5] += 0.5*v[5]; + } + } + } + + if (vflag_atom) { + if (newton_bond || i < nlocal) { + vatom[i][0] += 0.5*v[0]; + vatom[i][1] += 0.5*v[1]; + vatom[i][2] += 0.5*v[2]; + vatom[i][3] += 0.5*v[3]; + vatom[i][4] += 0.5*v[4]; + vatom[i][5] += 0.5*v[5]; + } + if (newton_bond || j < nlocal) { + vatom[j][0] += 0.5*v[0]; + vatom[j][1] += 0.5*v[1]; + vatom[j][2] += 0.5*v[2]; + vatom[j][3] += 0.5*v[3]; + vatom[j][4] += 0.5*v[4]; + vatom[j][5] += 0.5*v[5]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +double Bond::memory_usage() +{ + double bytes = maxeatom * sizeof(double); + bytes += maxvatom*6 * sizeof(double); + return bytes; +} diff --git a/src/dihedral.cpp b/src/dihedral.cpp new file mode 100644 index 0000000000..8f8c700e08 --- /dev/null +++ b/src/dihedral.cpp @@ -0,0 +1,249 @@ +/* ---------------------------------------------------------------------- + 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 "math.h" +#include "dihedral.h" +#include "atom.h" +#include "atom.h" +#include "force.h" +#include "pair.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + set dihedral contribution to Vdwl and Coulombic energy to 0.0 + DihedralCharmm will override this +------------------------------------------------------------------------- */ + +Dihedral::Dihedral(LAMMPS *lmp) : Pointers(lmp) +{ + energy = 0.0; + + allocated = 0; + PI = 4.0*atan(1.0); + + maxeatom = maxvatom = 0; + eatom = NULL; + vatom = NULL; +} + +/* ---------------------------------------------------------------------- */ + +Dihedral::~Dihedral() +{ + memory->sfree(eatom); + memory->destroy_2d_double_array(vatom); +} + +/* ---------------------------------------------------------------------- + check if all coeffs are set +------------------------------------------------------------------------- */ + +void Dihedral::init() +{ + if (!allocated) error->all("Dihedral coeffs are not set"); + for (int i = 1; i <= atom->ndihedraltypes; i++) + if (setflag[i] == 0) error->all("All dihedral coeffs are not set"); + init_style(); +} + +/* ---------------------------------------------------------------------- + setup for energy, virial computation + see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) +------------------------------------------------------------------------- */ + +void Dihedral::ev_setup(int eflag, int vflag) +{ + int i,n; + + evflag = 1; + + eflag_either = eflag; + eflag_global = eflag % 2; + eflag_atom = eflag / 2; + + vflag_either = vflag; + vflag_global = vflag % 4; + vflag_atom = vflag / 4; + + // reallocate per-atom arrays if necessary + + if (eflag_atom && atom->nmax > maxeatom) { + maxeatom = atom->nmax; + memory->sfree(eatom); + eatom = (double *) memory->smalloc(maxeatom*sizeof(double),"bond:eatom"); + } + if (vflag_atom && atom->nmax > maxvatom) { + maxvatom = atom->nmax; + memory->destroy_2d_double_array(vatom); + vatom = memory->create_2d_double_array(maxvatom,6,"bond:vatom"); + } + + // zero accumulators + + if (eflag_global) energy = 0.0; + if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; + if (eflag_atom) { + n = atom->nlocal; + if (force->newton_bond) n += atom->nghost; + for (i = 0; i < n; i++) eatom[i] = 0.0; + } + if (vflag_atom) { + n = atom->nlocal; + if (force->newton_bond) n += atom->nghost; + for (i = 0; i < n; i++) { + vatom[i][0] = 0.0; + vatom[i][1] = 0.0; + vatom[i][2] = 0.0; + vatom[i][3] = 0.0; + vatom[i][4] = 0.0; + vatom[i][5] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + tally energy and virial into global and per-atom accumulators + virial = r1F1 + r2F2 + r3F3 + r4F4 = (r1-r2) F1 + (r3-r2) F3 + (r4-r2) F4 + = (r1-r2) F1 + (r3-r2) F3 + (r4-r3 + r3-r2) F4 + = vb1*f1 + vb2*f3 + (vb3+vb2)*f4 +------------------------------------------------------------------------- */ + +void Dihedral::ev_tally(int i1, int i2, int i3, int i4, + int nlocal, int newton_bond, + double edihedral, double *f1, double *f3, double *f4, + double vb1x, double vb1y, double vb1z, + double vb2x, double vb2y, double vb2z, + double vb3x, double vb3y, double vb3z) +{ + double edihedralquarter,v[6]; + + if (eflag_either) { + if (eflag_global) { + if (newton_bond) energy += edihedral; + else { + edihedralquarter = 0.25*edihedral; + if (i1 < nlocal) energy += edihedralquarter; + if (i2 < nlocal) energy += edihedralquarter; + if (i3 < nlocal) energy += edihedralquarter; + if (i4 < nlocal) energy += edihedralquarter; + } + } + if (eflag_atom) { + edihedralquarter = 0.25*edihedral; + if (newton_bond || i1 < nlocal) eatom[i1] += edihedralquarter; + if (newton_bond || i2 < nlocal) eatom[i2] += edihedralquarter; + if (newton_bond || i3 < nlocal) eatom[i3] += edihedralquarter; + if (newton_bond || i4 < nlocal) eatom[i4] += edihedralquarter; + } + } + + if (vflag_either) { + v[0] = vb1x*f1[0] + vb2x*f3[0] + (vb3x+vb2x)*f4[0]; + v[1] = vb1y*f1[1] + vb2y*f3[1] + (vb3y+vb2y)*f4[1]; + v[2] = vb1z*f1[2] + vb2z*f3[2] + (vb3z+vb2z)*f4[2]; + v[3] = vb1x*f1[1] + vb2x*f3[1] + (vb3x+vb2x)*f4[1]; + v[4] = vb1x*f1[2] + vb2x*f3[2] + (vb3x+vb2x)*f4[2]; + v[5] = vb1y*f1[2] + vb2y*f3[2] + (vb3y+vb2y)*f4[2]; + + if (vflag_global) { + if (newton_bond) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } else { + if (i1 < nlocal) { + virial[0] += 0.25*v[0]; + virial[1] += 0.25*v[1]; + virial[2] += 0.25*v[2]; + virial[3] += 0.25*v[3]; + virial[4] += 0.25*v[4]; + virial[5] += 0.25*v[5]; + } + if (i2 < nlocal) { + virial[0] += 0.25*v[0]; + virial[1] += 0.25*v[1]; + virial[2] += 0.25*v[2]; + virial[3] += 0.25*v[3]; + virial[4] += 0.25*v[4]; + virial[5] += 0.25*v[5]; + } + if (i3 < nlocal) { + virial[0] += 0.25*v[0]; + virial[1] += 0.25*v[1]; + virial[2] += 0.25*v[2]; + virial[3] += 0.25*v[3]; + virial[4] += 0.25*v[4]; + virial[5] += 0.25*v[5]; + } + if (i4 < nlocal) { + virial[0] += 0.25*v[0]; + virial[1] += 0.25*v[1]; + virial[2] += 0.25*v[2]; + virial[3] += 0.25*v[3]; + virial[4] += 0.25*v[4]; + virial[5] += 0.25*v[5]; + } + } + } + + if (vflag_atom) { + if (newton_bond || i1 < nlocal) { + vatom[i1][0] += 0.25*v[0]; + vatom[i1][1] += 0.25*v[1]; + vatom[i1][2] += 0.25*v[2]; + vatom[i1][3] += 0.25*v[3]; + vatom[i1][4] += 0.25*v[4]; + vatom[i1][5] += 0.25*v[5]; + } + if (newton_bond || i2 < nlocal) { + vatom[i2][0] += 0.25*v[0]; + vatom[i2][1] += 0.25*v[1]; + vatom[i2][2] += 0.25*v[2]; + vatom[i2][3] += 0.25*v[3]; + vatom[i2][4] += 0.25*v[4]; + vatom[i2][5] += 0.25*v[5]; + } + if (newton_bond || i3 < nlocal) { + vatom[i3][0] += 0.25*v[0]; + vatom[i3][1] += 0.25*v[1]; + vatom[i3][2] += 0.25*v[2]; + vatom[i3][3] += 0.25*v[3]; + vatom[i3][4] += 0.25*v[4]; + vatom[i3][5] += 0.25*v[5]; + } + if (newton_bond || i4 < nlocal) { + vatom[i4][0] += 0.25*v[0]; + vatom[i4][1] += 0.25*v[1]; + vatom[i4][2] += 0.25*v[2]; + vatom[i4][3] += 0.25*v[3]; + vatom[i4][4] += 0.25*v[4]; + vatom[i4][5] += 0.25*v[5]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +double Dihedral::memory_usage() +{ + double bytes = maxeatom * sizeof(double); + bytes += maxvatom*6 * sizeof(double); + return bytes; +} diff --git a/src/improper.cpp b/src/improper.cpp new file mode 100644 index 0000000000..8697d9f722 --- /dev/null +++ b/src/improper.cpp @@ -0,0 +1,243 @@ +/* ---------------------------------------------------------------------- + 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 "math.h" +#include "improper.h" +#include "atom.h" +#include "force.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +Improper::Improper(LAMMPS *lmp) : Pointers(lmp) +{ + energy = 0.0; + + allocated = 0; + PI = 4.0*atan(1.0); + + maxeatom = maxvatom = 0; + eatom = NULL; + vatom = NULL; +} + +/* ---------------------------------------------------------------------- */ + +Improper::~Improper() +{ + memory->sfree(eatom); + memory->destroy_2d_double_array(vatom); +} + +/* ---------------------------------------------------------------------- + check if all coeffs are set +------------------------------------------------------------------------- */ + +void Improper::init() +{ + if (!allocated) error->all("Improper coeffs are not set"); + for (int i = 1; i <= atom->nimpropertypes; i++) + if (setflag[i] == 0) error->all("All improper coeffs are not set"); +} + +/* ---------------------------------------------------------------------- + setup for energy, virial computation + see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) +------------------------------------------------------------------------- */ + +void Improper::ev_setup(int eflag, int vflag) +{ + int i,n; + + evflag = 1; + + eflag_either = eflag; + eflag_global = eflag % 2; + eflag_atom = eflag / 2; + + vflag_either = vflag; + vflag_global = vflag % 4; + vflag_atom = vflag / 4; + + // reallocate per-atom arrays if necessary + + if (eflag_atom && atom->nmax > maxeatom) { + maxeatom = atom->nmax; + memory->sfree(eatom); + eatom = (double *) memory->smalloc(maxeatom*sizeof(double),"bond:eatom"); + } + if (vflag_atom && atom->nmax > maxvatom) { + maxvatom = atom->nmax; + memory->destroy_2d_double_array(vatom); + vatom = memory->create_2d_double_array(maxvatom,6,"bond:vatom"); + } + + // zero accumulators + + if (eflag_global) energy = 0.0; + if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; + if (eflag_atom) { + n = atom->nlocal; + if (force->newton_bond) n += atom->nghost; + for (i = 0; i < n; i++) eatom[i] = 0.0; + } + if (vflag_atom) { + n = atom->nlocal; + if (force->newton_bond) n += atom->nghost; + for (i = 0; i < n; i++) { + vatom[i][0] = 0.0; + vatom[i][1] = 0.0; + vatom[i][2] = 0.0; + vatom[i][3] = 0.0; + vatom[i][4] = 0.0; + vatom[i][5] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + tally energy and virial into global and per-atom accumulators + virial = r1F1 + r2F2 + r3F3 + r4F4 = (r1-r2) F1 + (r3-r2) F3 + (r4-r2) F4 + = (r1-r2) F1 + (r3-r2) F3 + (r4-r3 + r3-r2) F4 + = vb1*f1 + vb2*f3 + (vb3+vb2)*f4 +------------------------------------------------------------------------- */ + +void Improper::ev_tally(int i1, int i2, int i3, int i4, + int nlocal, int newton_bond, + double eimproper, double *f1, double *f3, double *f4, + double vb1x, double vb1y, double vb1z, + double vb2x, double vb2y, double vb2z, + double vb3x, double vb3y, double vb3z) +{ + double eimproperquarter,v[6]; + + if (eflag_either) { + if (eflag_global) { + if (newton_bond) energy += eimproper; + else { + eimproperquarter = 0.25*eimproper; + if (i1 < nlocal) energy += eimproperquarter; + if (i2 < nlocal) energy += eimproperquarter; + if (i3 < nlocal) energy += eimproperquarter; + if (i4 < nlocal) energy += eimproperquarter; + } + } + if (eflag_atom) { + eimproperquarter = 0.25*eimproper; + if (newton_bond || i1 < nlocal) eatom[i1] += eimproperquarter; + if (newton_bond || i2 < nlocal) eatom[i2] += eimproperquarter; + if (newton_bond || i3 < nlocal) eatom[i3] += eimproperquarter; + if (newton_bond || i4 < nlocal) eatom[i4] += eimproperquarter; + } + } + + if (vflag_either) { + v[0] = vb1x*f1[0] + vb2x*f3[0] + (vb3x+vb2x)*f4[0]; + v[1] = vb1y*f1[1] + vb2y*f3[1] + (vb3y+vb2y)*f4[1]; + v[2] = vb1z*f1[2] + vb2z*f3[2] + (vb3z+vb2z)*f4[2]; + v[3] = vb1x*f1[1] + vb2x*f3[1] + (vb3x+vb2x)*f4[1]; + v[4] = vb1x*f1[2] + vb2x*f3[2] + (vb3x+vb2x)*f4[2]; + v[5] = vb1y*f1[2] + vb2y*f3[2] + (vb3y+vb2y)*f4[2]; + + if (vflag_global) { + if (newton_bond) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } else { + if (i1 < nlocal) { + virial[0] += 0.25*v[0]; + virial[1] += 0.25*v[1]; + virial[2] += 0.25*v[2]; + virial[3] += 0.25*v[3]; + virial[4] += 0.25*v[4]; + virial[5] += 0.25*v[5]; + } + if (i2 < nlocal) { + virial[0] += 0.25*v[0]; + virial[1] += 0.25*v[1]; + virial[2] += 0.25*v[2]; + virial[3] += 0.25*v[3]; + virial[4] += 0.25*v[4]; + virial[5] += 0.25*v[5]; + } + if (i3 < nlocal) { + virial[0] += 0.25*v[0]; + virial[1] += 0.25*v[1]; + virial[2] += 0.25*v[2]; + virial[3] += 0.25*v[3]; + virial[4] += 0.25*v[4]; + virial[5] += 0.25*v[5]; + } + if (i4 < nlocal) { + virial[0] += 0.25*v[0]; + virial[1] += 0.25*v[1]; + virial[2] += 0.25*v[2]; + virial[3] += 0.25*v[3]; + virial[4] += 0.25*v[4]; + virial[5] += 0.25*v[5]; + } + } + } + + if (vflag_atom) { + if (newton_bond || i1 < nlocal) { + vatom[i1][0] += 0.25*v[0]; + vatom[i1][1] += 0.25*v[1]; + vatom[i1][2] += 0.25*v[2]; + vatom[i1][3] += 0.25*v[3]; + vatom[i1][4] += 0.25*v[4]; + vatom[i1][5] += 0.25*v[5]; + } + if (newton_bond || i2 < nlocal) { + vatom[i2][0] += 0.25*v[0]; + vatom[i2][1] += 0.25*v[1]; + vatom[i2][2] += 0.25*v[2]; + vatom[i2][3] += 0.25*v[3]; + vatom[i2][4] += 0.25*v[4]; + vatom[i2][5] += 0.25*v[5]; + } + if (newton_bond || i3 < nlocal) { + vatom[i3][0] += 0.25*v[0]; + vatom[i3][1] += 0.25*v[1]; + vatom[i3][2] += 0.25*v[2]; + vatom[i3][3] += 0.25*v[3]; + vatom[i3][4] += 0.25*v[4]; + vatom[i3][5] += 0.25*v[5]; + } + if (newton_bond || i4 < nlocal) { + vatom[i4][0] += 0.25*v[0]; + vatom[i4][1] += 0.25*v[1]; + vatom[i4][2] += 0.25*v[2]; + vatom[i4][3] += 0.25*v[3]; + vatom[i4][4] += 0.25*v[4]; + vatom[i4][5] += 0.25*v[5]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +double Improper::memory_usage() +{ + double bytes = maxeatom * sizeof(double); + bytes += maxvatom*6 * sizeof(double); + return bytes; +}