From 94ddc19477beb8379e74f8db745ba6b63b2763f4 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 12 Feb 2009 21:36:52 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2589 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/COLLOID/pair_colloid.cpp | 8 - src/COLLOID/pair_colloid.h | 1 - src/GRANULAR/atom_vec_granular.cpp | 2 + src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp | 12 +- src/MANYBODY/pair_sw.cpp | 5 +- src/MANYBODY/pair_tersoff.cpp | 5 +- src/PERI/atom_vec_peri.cpp | 1220 ++++++++++---------- src/PERI/pair_peri_pmb.cpp | 1042 ++++++++--------- src/atom.h | 1 + src/compute_centro_atom.h | 2 +- src/compute_coord_atom.cpp | 9 +- src/compute_coord_atom.h | 4 +- src/compute_erotate_sphere.cpp | 16 +- src/compute_ke.cpp | 11 +- src/compute_ke_atom.cpp | 6 +- src/compute_stress_atom.cpp | 6 +- src/compute_temp.cpp | 14 +- src/compute_temp_com.cpp | 12 +- src/compute_temp_deform.cpp | 15 +- src/compute_temp_partial.cpp | 16 +- src/compute_temp_ramp.cpp | 12 +- src/compute_temp_region.cpp | 12 +- src/compute_temp_sphere.cpp | 56 +- src/dump_custom.cpp | 14 +- src/fix_ave_spatial.cpp | 4 +- src/fix_dt_reset.cpp | 4 +- src/fix_gravity.cpp | 11 +- src/fix_line_force.cpp | 2 - src/fix_nve.cpp | 12 +- src/fix_nve_limit.cpp | 12 +- src/fix_nve_sphere.cpp | 76 +- src/fix_shake.cpp | 213 +++- src/fix_shake.h | 5 +- src/fix_thermal_conductivity.cpp | 8 +- src/fix_viscosity.cpp | 19 +- src/group.cpp | 44 +- src/pair_lj_cut.cpp | 8 - src/pair_lj_cut.h | 1 - src/velocity.cpp | 12 +- 39 files changed, 1506 insertions(+), 1426 deletions(-) diff --git a/src/COLLOID/pair_colloid.cpp b/src/COLLOID/pair_colloid.cpp index a98f1fa8f9..8d119b046b 100644 --- a/src/COLLOID/pair_colloid.cpp +++ b/src/COLLOID/pair_colloid.cpp @@ -521,11 +521,3 @@ double PairColloid::single(int i, int j, int itype, int jtype, double rsq, return factor_lj*phi; } - -/* ---------------------------------------------------------------------- */ - -void *PairColloid::extract(char *str) -{ - if (strcmp(str,"diameter") == 0) return (void *) diameter; - return NULL; -} diff --git a/src/COLLOID/pair_colloid.h b/src/COLLOID/pair_colloid.h index 00f41b4b36..976f1e6f92 100644 --- a/src/COLLOID/pair_colloid.h +++ b/src/COLLOID/pair_colloid.h @@ -31,7 +31,6 @@ class PairColloid : public Pair { void write_restart_settings(FILE *); void read_restart_settings(FILE *); double single(int, int, int, int, double, double, double, double &); - void *extract(char *); private: double cut_global; diff --git a/src/GRANULAR/atom_vec_granular.cpp b/src/GRANULAR/atom_vec_granular.cpp index 9d313e8a47..b77513a64a 100644 --- a/src/GRANULAR/atom_vec_granular.cpp +++ b/src/GRANULAR/atom_vec_granular.cpp @@ -604,6 +604,7 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values) density[nlocal] = atof(values[3]); rmass[nlocal] = 4.0*PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; + if (rmass[nlocal] <= 0.0) error->one("Invalid mass value"); x[nlocal][0] = coord[0]; x[nlocal][1] = coord[1]; @@ -633,6 +634,7 @@ int AtomVecGranular::data_atom_hybrid(int nlocal, char **values) density[nlocal] = atof(values[1]); rmass[nlocal] = 4.0*PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; + if (rmass[nlocal] <= 0.0) error->one("Invalid mass value"); v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; diff --git a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp index fa8c664448..abf1798f6d 100644 --- a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp @@ -54,6 +54,10 @@ PairLJCutCoulLongTIP4P::PairLJCutCoulLongTIP4P(LAMMPS *lmp) : PairLJCutCoulLong(lmp) { single_enable = 0; + + // TIP4P cannot compute virial as F dot r + // due to find_M() finding bonded H atoms which are not near O atom + no_virial_compute = 1; } @@ -76,15 +80,9 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) float rsq; int *int_rsq = (int *) &rsq; - // if global component of incoming vflag = 2, reset vflag as if it were 1 - // necessary since TIP4P cannot compute virial as F dot r - // due to find_M() finding bonded H atoms which are not near O atom - - if (vflag % 4 == 2) vflag = 1 + vflag/4 * 4; - evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = 0; + else evflag = vflag_fdotr = 0; double **f = atom->f; double **x = atom->x; diff --git a/src/MANYBODY/pair_sw.cpp b/src/MANYBODY/pair_sw.cpp index c590021225..886f7abf52 100755 --- a/src/MANYBODY/pair_sw.cpp +++ b/src/MANYBODY/pair_sw.cpp @@ -122,9 +122,8 @@ void PairSW::compute(int eflag, int vflag) if ((itag+jtag) % 2 == 1) continue; } else { if (x[j][2] < ztmp) continue; - else if (x[j][2] == ztmp && x[j][1] < ytmp) continue; - else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) - continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; } jtype = map[type[j]]; diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index 3a70e6cd63..cfa133d3d7 100755 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -125,9 +125,8 @@ void PairTersoff::compute(int eflag, int vflag) if ((itag+jtag) % 2 == 1) continue; } else { if (x[j][2] < x[i][2]) continue; - else if (x[j][2] == ztmp && x[j][1] < ytmp) continue; - else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) - continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; } jtype = map[type[j]]; diff --git a/src/PERI/atom_vec_peri.cpp b/src/PERI/atom_vec_peri.cpp index 9c23b75721..4fdf661f5b 100644 --- a/src/PERI/atom_vec_peri.cpp +++ b/src/PERI/atom_vec_peri.cpp @@ -1,610 +1,610 @@ -/* ---------------------------------------------------------------------- - 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: Mike Parks (SNL) -------------------------------------------------------------------------- */ - -#include "float.h" -#include "stdlib.h" -#include "atom_vec_peri.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 - -/* ---------------------------------------------------------------------- */ - -AtomVecPeri::AtomVecPeri(LAMMPS *lmp, int narg, char **arg) : -AtomVec(lmp, narg, arg) -{ - comm_x_only = 0; - size_comm = 4; - size_reverse = 3; - size_border = 11; - size_data_atom = 7; - size_data_vel = 4; - xcol_data = 5; - - atom->vfrac_flag = atom->density_flag = atom->rmass_flag = 1; -} - -/* ---------------------------------------------------------------------- - grow atom arrays - n = 0 grows arrays by DELTA - n > 0 allocates arrays to size n -------------------------------------------------------------------------- */ - -void AtomVecPeri::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"); - - vfrac = atom->vfrac = (double *) - memory->srealloc(atom->vfrac,nmax*sizeof(double),"atom:vfrac"); - density = atom->density = (double *) - memory->srealloc(atom->density,nmax*sizeof(double),"atom:density"); - rmass = atom->rmass = (double *) - memory->srealloc(atom->rmass,nmax*sizeof(double),"atom:rmass"); - s0 = atom->s0 = (double *) - memory->srealloc(atom->s0,nmax*sizeof(double),"atom:s0"); - x0 = atom->x0 = memory->grow_2d_double_array(atom->x0,nmax,3,"atom:x0"); - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecPeri::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]; - - vfrac[j] = vfrac[i]; - density[j] = density[i]; - rmass[j] = rmass[i]; - s0[j] = s0[i]; - x0[j][0] = x0[i][0]; - x0[j][1] = x0[i][1]; - x0[j][2] = x0[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 AtomVecPeri::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++] = s0[j]; - } - } 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++] = s0[j]; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecPeri::pack_comm_one(int i, double *buf) -{ - buf[0] = s0[i]; - return 1; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecPeri::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++]; - s0[i] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecPeri::unpack_comm_one(int i, double *buf) -{ - s0[i] = buf[0]; - return 1; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecPeri::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]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecPeri::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++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecPeri::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++] = vfrac[j]; - buf[m++] = s0[j]; - buf[m++] = x0[j][0]; - buf[m++] = x0[j][1]; - buf[m++] = x0[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++] = vfrac[j]; - buf[m++] = s0[j]; - buf[m++] = x0[j][0]; - buf[m++] = x0[j][1]; - buf[m++] = x0[j][2]; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecPeri::pack_border_one(int i, double *buf) -{ - buf[0] = vfrac[i]; - buf[1] = s0[i]; - buf[2] = x0[i][0]; - buf[3] = x0[i][1]; - buf[4] = x0[i][2]; - return 5; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecPeri::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++]); - vfrac[i] = buf[m++]; - s0[i] = buf[m++]; - x0[i][0] = buf[m++]; - x0[i][1] = buf[m++]; - x0[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecPeri::unpack_border_one(int i, double *buf) -{ - vfrac[i] = buf[0]; - s0[i] = buf[1]; - x0[i][0] = buf[2]; - x0[i][1] = buf[3]; - x0[i][2] = buf[4]; - return 5; -} - -/* ---------------------------------------------------------------------- - pack data for atom I for sending to another proc - xyz must be 1st 3 values, so comm::exchange() can test on them -------------------------------------------------------------------------- */ - -int AtomVecPeri::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++] = vfrac[i]; - buf[m++] = density[i]; - buf[m++] = rmass[i]; - buf[m++] = s0[i]; - buf[m++] = x0[i][0]; - buf[m++] = x0[i][1]; - buf[m++] = x0[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 AtomVecPeri::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++]); - - vfrac[nlocal] = buf[m++]; - density[nlocal] = buf[m++]; - rmass[nlocal] = buf[m++]; - s0[nlocal] = buf[m++]; - x0[nlocal][0] = buf[m++]; - x0[nlocal][1] = buf[m++]; - x0[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 AtomVecPeri::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 AtomVecPeri::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++] = vfrac[i]; - buf[m++] = density[i]; - buf[m++] = rmass[i]; - buf[m++] = s0[i]; - buf[m++] = x0[i][0]; - buf[m++] = x0[i][1]; - buf[m++] = x0[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 AtomVecPeri::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++]; - - vfrac[nlocal] = buf[m++]; - density[nlocal] = buf[m++]; - rmass[nlocal] = buf[m++]; - s0[nlocal] = buf[m++]; - x0[nlocal][0] = buf[m++]; - x0[nlocal][1] = buf[m++]; - x0[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 AtomVecPeri::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; - - vfrac[nlocal] = 1.0; - density[nlocal] = 1.0; - rmass[nlocal] = density[nlocal]; - s0[nlocal] = DBL_MAX; - x0[nlocal][0] = coord[0]; - x0[nlocal][1] = coord[1]; - x0[nlocal][2] = coord[2]; - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - unpack one line from Atoms section of data file - initialize other atom quantities -------------------------------------------------------------------------- */ - -void AtomVecPeri::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"); - - vfrac[nlocal] = atof(values[2]); - density[nlocal] = atof(values[3]); - rmass[nlocal] = density[nlocal]; - - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - - image[nlocal] = imagetmp; - - mask[nlocal] = 1; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - - s0[nlocal] = DBL_MAX; - x0[nlocal][0] = coord[0]; - x0[nlocal][1] = coord[1]; - x0[nlocal][2] = coord[2]; - - atom->nlocal++; -} - - -/* ---------------------------------------------------------------------- - unpack hybrid quantities from one line in Atoms section of data file - initialize other atom quantities for this sub-style -------------------------------------------------------------------------- */ - -int AtomVecPeri::data_atom_hybrid(int nlocal, char **values) -{ - vfrac[nlocal] = atof(values[0]); - density[nlocal] = atof(values[1]); - rmass[nlocal] = density[nlocal]; - - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - - s0[nlocal] = DBL_MAX; - x0[nlocal][0] = x[nlocal][0]; - x0[nlocal][1] = x[nlocal][1]; - x0[nlocal][2] = x[nlocal][2]; - - return 2; -} - -/* ---------------------------------------------------------------------- - return # of bytes of allocated memory -------------------------------------------------------------------------- */ - -double AtomVecPeri::memory_usage() -{ - double bytes = 0.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("vfrac")) bytes += nmax * sizeof(double); - if (atom->memcheck("density")) bytes += nmax * sizeof(double); - if (atom->memcheck("rmass")) bytes += nmax * sizeof(double); - if (atom->memcheck("s0")) bytes += nmax * sizeof(double); - if (atom->memcheck("x0")) bytes += nmax*3 * sizeof(double); - - return bytes; -} +/* ---------------------------------------------------------------------- + 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: Mike Parks (SNL) +------------------------------------------------------------------------- */ + +#include "float.h" +#include "stdlib.h" +#include "atom_vec_peri.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 + +/* ---------------------------------------------------------------------- */ + +AtomVecPeri::AtomVecPeri(LAMMPS *lmp, int narg, char **arg) : +AtomVec(lmp, narg, arg) +{ + comm_x_only = 0; + size_comm = 4; + size_reverse = 3; + size_border = 11; + size_data_atom = 7; + size_data_vel = 4; + xcol_data = 5; + + atom->vfrac_flag = atom->density_flag = atom->rmass_flag = 1; +} + +/* ---------------------------------------------------------------------- + grow atom arrays + n = 0 grows arrays by DELTA + n > 0 allocates arrays to size n +------------------------------------------------------------------------- */ + +void AtomVecPeri::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"); + + vfrac = atom->vfrac = (double *) + memory->srealloc(atom->vfrac,nmax*sizeof(double),"atom:vfrac"); + density = atom->density = (double *) + memory->srealloc(atom->density,nmax*sizeof(double),"atom:density"); + rmass = atom->rmass = (double *) + memory->srealloc(atom->rmass,nmax*sizeof(double),"atom:rmass"); + s0 = atom->s0 = (double *) + memory->srealloc(atom->s0,nmax*sizeof(double),"atom:s0"); + x0 = atom->x0 = memory->grow_2d_double_array(atom->x0,nmax,3,"atom:x0"); + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecPeri::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]; + + vfrac[j] = vfrac[i]; + density[j] = density[i]; + rmass[j] = rmass[i]; + s0[j] = s0[i]; + x0[j][0] = x0[i][0]; + x0[j][1] = x0[i][1]; + x0[j][2] = x0[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 AtomVecPeri::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++] = s0[j]; + } + } 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++] = s0[j]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::pack_comm_one(int i, double *buf) +{ + buf[0] = s0[i]; + return 1; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecPeri::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++]; + s0[i] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::unpack_comm_one(int i, double *buf) +{ + s0[i] = buf[0]; + return 1; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::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]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecPeri::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++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::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++] = vfrac[j]; + buf[m++] = s0[j]; + buf[m++] = x0[j][0]; + buf[m++] = x0[j][1]; + buf[m++] = x0[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++] = vfrac[j]; + buf[m++] = s0[j]; + buf[m++] = x0[j][0]; + buf[m++] = x0[j][1]; + buf[m++] = x0[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::pack_border_one(int i, double *buf) +{ + buf[0] = vfrac[i]; + buf[1] = s0[i]; + buf[2] = x0[i][0]; + buf[3] = x0[i][1]; + buf[4] = x0[i][2]; + return 5; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecPeri::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++]); + vfrac[i] = buf[m++]; + s0[i] = buf[m++]; + x0[i][0] = buf[m++]; + x0[i][1] = buf[m++]; + x0[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::unpack_border_one(int i, double *buf) +{ + vfrac[i] = buf[0]; + s0[i] = buf[1]; + x0[i][0] = buf[2]; + x0[i][1] = buf[3]; + x0[i][2] = buf[4]; + return 5; +} + +/* ---------------------------------------------------------------------- + pack data for atom I for sending to another proc + xyz must be 1st 3 values, so comm::exchange() can test on them +------------------------------------------------------------------------- */ + +int AtomVecPeri::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++] = vfrac[i]; + buf[m++] = density[i]; + buf[m++] = rmass[i]; + buf[m++] = s0[i]; + buf[m++] = x0[i][0]; + buf[m++] = x0[i][1]; + buf[m++] = x0[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 AtomVecPeri::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++]); + + vfrac[nlocal] = buf[m++]; + density[nlocal] = buf[m++]; + rmass[nlocal] = buf[m++]; + s0[nlocal] = buf[m++]; + x0[nlocal][0] = buf[m++]; + x0[nlocal][1] = buf[m++]; + x0[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 AtomVecPeri::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 AtomVecPeri::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++] = vfrac[i]; + buf[m++] = density[i]; + buf[m++] = rmass[i]; + buf[m++] = s0[i]; + buf[m++] = x0[i][0]; + buf[m++] = x0[i][1]; + buf[m++] = x0[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 AtomVecPeri::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++]; + + vfrac[nlocal] = buf[m++]; + density[nlocal] = buf[m++]; + rmass[nlocal] = buf[m++]; + s0[nlocal] = buf[m++]; + x0[nlocal][0] = buf[m++]; + x0[nlocal][1] = buf[m++]; + x0[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 AtomVecPeri::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; + + vfrac[nlocal] = 1.0; + density[nlocal] = 1.0; + rmass[nlocal] = density[nlocal]; + s0[nlocal] = DBL_MAX; + x0[nlocal][0] = coord[0]; + x0[nlocal][1] = coord[1]; + x0[nlocal][2] = coord[2]; + + atom->nlocal++; +} + +/* ---------------------------------------------------------------------- + unpack one line from Atoms section of data file + initialize other atom quantities +------------------------------------------------------------------------- */ + +void AtomVecPeri::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"); + + vfrac[nlocal] = atof(values[2]); + density[nlocal] = atof(values[3]); + rmass[nlocal] = density[nlocal]; + + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + + image[nlocal] = imagetmp; + + mask[nlocal] = 1; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + + s0[nlocal] = DBL_MAX; + x0[nlocal][0] = coord[0]; + x0[nlocal][1] = coord[1]; + x0[nlocal][2] = coord[2]; + + atom->nlocal++; +} + + +/* ---------------------------------------------------------------------- + unpack hybrid quantities from one line in Atoms section of data file + initialize other atom quantities for this sub-style +------------------------------------------------------------------------- */ + +int AtomVecPeri::data_atom_hybrid(int nlocal, char **values) +{ + vfrac[nlocal] = atof(values[0]); + density[nlocal] = atof(values[1]); + rmass[nlocal] = density[nlocal]; + + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + + s0[nlocal] = DBL_MAX; + x0[nlocal][0] = x[nlocal][0]; + x0[nlocal][1] = x[nlocal][1]; + x0[nlocal][2] = x[nlocal][2]; + + return 2; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory +------------------------------------------------------------------------- */ + +double AtomVecPeri::memory_usage() +{ + double bytes = 0.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("vfrac")) bytes += nmax * sizeof(double); + if (atom->memcheck("density")) bytes += nmax * sizeof(double); + if (atom->memcheck("rmass")) bytes += nmax * sizeof(double); + if (atom->memcheck("s0")) bytes += nmax * sizeof(double); + if (atom->memcheck("x0")) bytes += nmax*3 * sizeof(double); + + return bytes; +} diff --git a/src/PERI/pair_peri_pmb.cpp b/src/PERI/pair_peri_pmb.cpp index c21c4aaad2..8a801ba27e 100644 --- a/src/PERI/pair_peri_pmb.cpp +++ b/src/PERI/pair_peri_pmb.cpp @@ -1,521 +1,521 @@ -/* ---------------------------------------------------------------------- - 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: Mike Parks (SNL) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "float.h" -#include "stdlib.h" -#include "string.h" -#include "pair_peri_pmb.h" -#include "atom.h" -#include "domain.h" -#include "lattice.h" -#include "force.h" -#include "update.h" -#include "modify.h" -#include "fix.h" -#include "fix_peri_neigh.h" -#include "comm.h" -#include "neighbor.h" -#include "neigh_list.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)) - -/* ---------------------------------------------------------------------- */ - -PairPeriPMB::PairPeriPMB(LAMMPS *lmp) : Pair(lmp) -{ - for (int i = 0; i < 6; i++) virial[i] = 0.0; - - ifix_peri = -1; - - nmax = 0; - s0_new = NULL; - - kspring = NULL; - s00 = NULL; - alpha = NULL; - cut = NULL; -} - -/* ---------------------------------------------------------------------- */ - -PairPeriPMB::~PairPeriPMB() -{ - if (ifix_peri >= 0) modify->delete_fix("PERI_NEIGH"); - - if (allocated) { - memory->destroy_2d_int_array(setflag); - memory->destroy_2d_double_array(cutsq); - memory->destroy_2d_double_array(kspring); - memory->destroy_2d_double_array(s00); - memory->destroy_2d_double_array(alpha); - memory->destroy_2d_double_array(cut); - memory->sfree(s0_new); - } -} - -/* ---------------------------------------------------------------------- */ - -void PairPeriPMB::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz; - double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0,rsq0; - double rsq,r,dr,rk,evdwl,fpair,fbond; - int *ilist,*jlist,*numneigh,**firstneigh; - double d_ij,delta,stretch; - - evdwl = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - double **f = atom->f; - double **x = atom->x; - int *type = atom->type; - int nlocal = atom->nlocal; - - double *vfrac = atom->vfrac; - double *s0 = atom->s0; - double **x0 = atom->x0; - double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; - int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; - int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; - - // lc = lattice constant - // init_style guarantees it's the same in x, y, and z - - double lc = domain->lattice->xlattice; - double half_lc = 0.5*lc; - double vfrac_scale = 1.0; - - // short-range forces - - int nall = atom->nlocal + atom->nghost; - int newton_pair = force->newton_pair; - int nonperiodic = domain->nonperiodic; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // loop over neighbors of my atoms - // need minimg() for x0 difference since not ghosted - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - xtmp0 = x0[i][0]; - ytmp0 = x0[i][1]; - ztmp0 = x0[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j %= nall; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - delx0 = xtmp0 - x0[j][0]; - dely0 = ytmp0 - x0[j][1]; - delz0 = ztmp0 - x0[j][2]; - if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0); - rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; - jtype = type[j]; - - r = sqrt(rsq); - - // short-range interaction distance based on initial particle position - // 0.9 and 1.35 are constants - - d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); - - // short-range contact forces - // 15 is constant taken from the EMU Theory Manual - // Silling, 12 May 2005, p 18 - - if (r < d_ij) { - dr = r - d_ij; - - rk = (15.0 * kspring[itype][jtype] * vfrac[j]) * - (dr / sqrt(cutsq[itype][jtype])); - if (r > 0.0) fpair = -(rk/r); - else fpair = 0.0; - - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; - } - - if (eflag) evdwl = 0.5*rk*dr; - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); - } - } - } - - // bond forces - - if (atom->nmax > nmax) { - memory->sfree(s0_new); - nmax = atom->nmax; - s0_new = (double *) memory->smalloc(nmax*sizeof(double),"pair:s0_new"); - } - - // loop over my particles and their partners - // partner list contains all bond partners, so I-J appears twice - // if bond already broken, skip this partner - // first = true if this is first neighbor of particle i - - bool first; - - for (i = 0; i < nlocal; i++) { - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jnum = npartner[i]; - s0_new[i] = DBL_MAX; - first = true; - - for (jj = 0; jj < jnum; jj++) { - if (partner[i][jj] == 0) continue; - j = atom->map(partner[i][jj]); - - // check if lost a partner without first breaking bond - - if (j < 0) { - partner[i][jj] = 0; - continue; - } - - // compute force density, add to PD equation of motion - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - if (nonperiodic == 0) domain->minimum_image(delx,dely,delz); - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - delta = sqrt(cutsq[itype][jtype]); - r = sqrt(rsq); - dr = r - r0[i][jj]; - - // avoid roundoff errors - - if (fabs(dr) < 2.2204e-016) dr = 0.0; - - // scale vfrac[j] if particle j near the horizon - - if ((fabs(r0[i][jj] - delta)) <= half_lc) - vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + - (1.0 + ((delta - half_lc)/(2*half_lc) ) ); - else vfrac_scale = 1.0; - - stretch = dr / r0[i][jj]; - rk = (kspring[itype][jtype] * vfrac[j]) * vfrac_scale * stretch; - if (r > 0.0) fbond = -(rk/r); - else fbond = 0.0; - - f[i][0] += delx*fbond; - f[i][1] += dely*fbond; - f[i][2] += delz*fbond; - - // since I-J is double counted, set newton off & use 1/2 factor and I,I - - if (eflag) evdwl = 0.5*rk*dr; - if (evflag) ev_tally(i,i,nlocal,0, - 0.5*evdwl,0.0,0.5*fbond,delx,dely,delz); - - // find stretch in bond I-J and break if necessary - // use s0 from previous timestep - - if (stretch > MIN(s0[i],s0[j])) partner[i][jj] = 0; - - // update s0 for next timestep - - if (first) - s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch); - else - s0_new[i] = MAX(s0_new[i], - s00[itype][jtype] - (alpha[itype][jtype] * stretch)); - first = false; - } - } - - if (vflag_fdotr) virial_compute(); - - // store new s0 - - for (i = 0; i < nlocal; i++) s0[i] = s0_new[i]; -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairPeriPMB::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"); - kspring = memory->create_2d_double_array(n+1,n+1,"pair:kspring"); - s00 = memory->create_2d_double_array(n+1,n+1,"pair:s00"); - alpha = memory->create_2d_double_array(n+1,n+1,"pair:alpha"); - cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairPeriPMB::settings(int narg, char **arg) -{ - if (narg) error->all("Illegal pair_style command"); -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairPeriPMB::coeff(int narg, char **arg) -{ - if (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 kspring_one = atof(arg[2]); - double cut_one = atof(arg[3]); - double s00_one = atof(arg[4]); - double alpha_one = atof(arg[5]); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - kspring[i][j] = kspring_one; - s00[i][j] = s00_one; - alpha[i][j] = alpha_one; - cut[i][j] = cut_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 PairPeriPMB::init_one(int i, int j) -{ - if (setflag[i][j] == 0) error->all("All pair coeffs are not set"); - - cutsq[i][j] = cut[i][j] * cut[i][j]; - cutsq[j][i] = cutsq[i][j]; - - // set other j,i parameters - - kspring[j][i] = kspring[i][j]; - alpha[j][i] = alpha[i][j]; - s00[j][i] = s00[i][j]; - - return cut[i][j]; -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairPeriPMB::init_style() -{ - // error checks - - if (atom->map_style == 0) - error->all("Pair peri requires an atom map, see atom_modify"); - - if (atom->style_match("peri") == 0) - error->all("Pair style peri_pmb requires atom style peri"); - - if (domain->lattice == NULL) - error->all("Pair peri requires a lattice be defined"); - if (domain->lattice->xlattice != domain->lattice->ylattice || - domain->lattice->xlattice != domain->lattice->zlattice || - domain->lattice->ylattice != domain->lattice->zlattice) - error->all("Pair peri lattice is not identical in x, y, and z"); - - // if first init, create Fix needed for storing fixed neighbors - - if (ifix_peri == -1) { - char **fixarg = new char*[3]; - fixarg[0] = (char *) "PERI_NEIGH"; - fixarg[1] = (char *) "all"; - fixarg[2] = (char *) "PERI_NEIGH"; - modify->add_fix(3,fixarg); - delete [] fixarg; - } - - // find associated PERI_NEIGH fix that must exist - // could have changed locations in fix list since created - - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"PERI_NEIGH") == 0) ifix_peri = i; - - int irequest = neighbor->request(this); -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairPeriPMB::write_restart(FILE *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(&kspring[i][j],sizeof(double),1,fp); - fwrite(&s00[i][j],sizeof(double),1,fp); - fwrite(&alpha[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairPeriPMB::read_restart(FILE *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(&kspring[i][j],sizeof(double),1,fp); - fread(&s00[i][j],sizeof(double),1,fp); - fread(&alpha[i][j],sizeof(double),1,fp); - fread(&cut[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&kspring[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&s00[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&alpha[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); - } - } -} - -/* ---------------------------------------------------------------------- */ - -double PairPeriPMB::single(int i, int j, int itype, int jtype, double rsq, - double factor_coul, double factor_lj, - double &fforce) -{ - double delx0,dely0,delz0,rsq0; - double d_ij,r,dr,rk,vfrac_scale; - - double *vfrac = atom->vfrac; - double **x0 = atom->x0; - double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; - int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; - int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; - - double lc = domain->lattice->xlattice; - double half_lc = 0.5*lc; - - delx0 = x0[i][0] - x0[j][0]; - dely0 = x0[i][1] - x0[j][1]; - delz0 = x0[i][2] - x0[j][2]; - if (domain->nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0); - rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; - - d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); - r = sqrt(rsq); - - double energy = 0.0; - fforce = 0.0; - - if (r < d_ij) { - dr = r - d_ij; - rk = (15.0 * kspring[itype][jtype] * vfrac[j]) * - (dr / sqrt(cutsq[itype][jtype])); - if (r > 0.0) fforce += -(rk/r); - energy += 0.5*rk*dr; - } - - int jnum = npartner[i]; - for (int jj = 0; jj < jnum; jj++) { - if (partner[i][jj] == 0) continue; - if (j < 0) continue; - if (j == atom->map(partner[i][jj])) { - dr = r - r0[i][jj]; - if (fabs(dr) < 2.2204e-016) dr = 0.0; - if ( (fabs(r0[i][jj] - sqrt(cutsq[itype][jtype]))) <= half_lc) - vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + - (1.0 + ((sqrt(cutsq[itype][jtype]) - half_lc)/(2*half_lc))); - else vfrac_scale = 1.0; - rk = (kspring[itype][jtype] * vfrac[j] * vfrac_scale) * - (dr / r0[i][jj]); - if (r > 0.0) fforce += -(rk/r); - energy += 0.5*rk*dr; - } - } - - return energy; -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based arrays -------------------------------------------------------------------------- */ - -double PairPeriPMB::memory_usage() -{ - double bytes = nmax * sizeof(double); - return bytes; -} +/* ---------------------------------------------------------------------- + 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: Mike Parks (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "float.h" +#include "stdlib.h" +#include "string.h" +#include "pair_peri_pmb.h" +#include "atom.h" +#include "domain.h" +#include "lattice.h" +#include "force.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "fix_peri_neigh.h" +#include "comm.h" +#include "neighbor.h" +#include "neigh_list.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)) + +/* ---------------------------------------------------------------------- */ + +PairPeriPMB::PairPeriPMB(LAMMPS *lmp) : Pair(lmp) +{ + for (int i = 0; i < 6; i++) virial[i] = 0.0; + + ifix_peri = -1; + + nmax = 0; + s0_new = NULL; + + kspring = NULL; + s00 = NULL; + alpha = NULL; + cut = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairPeriPMB::~PairPeriPMB() +{ + if (ifix_peri >= 0) modify->delete_fix("PERI_NEIGH"); + + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + memory->destroy_2d_double_array(kspring); + memory->destroy_2d_double_array(s00); + memory->destroy_2d_double_array(alpha); + memory->destroy_2d_double_array(cut); + memory->sfree(s0_new); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriPMB::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz; + double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0,rsq0; + double rsq,r,dr,rk,evdwl,fpair,fbond; + int *ilist,*jlist,*numneigh,**firstneigh; + double d_ij,delta,stretch; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **f = atom->f; + double **x = atom->x; + int *type = atom->type; + int nlocal = atom->nlocal; + + double *vfrac = atom->vfrac; + double *s0 = atom->s0; + double **x0 = atom->x0; + double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; + int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + + // lc = lattice constant + // init_style guarantees it's the same in x, y, and z + + double lc = domain->lattice->xlattice; + double half_lc = 0.5*lc; + double vfrac_scale = 1.0; + + // short-range forces + + int nall = atom->nlocal + atom->nghost; + int newton_pair = force->newton_pair; + int nonperiodic = domain->nonperiodic; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + // need minimg() for x0 difference since not ghosted + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j %= nall; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0); + rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; + jtype = type[j]; + + r = sqrt(rsq); + + // short-range interaction distance based on initial particle position + // 0.9 and 1.35 are constants + + d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); + + // short-range contact forces + // 15 is constant taken from the EMU Theory Manual + // Silling, 12 May 2005, p 18 + + if (r < d_ij) { + dr = r - d_ij; + + rk = (15.0 * kspring[itype][jtype] * vfrac[j]) * + (dr / sqrt(cutsq[itype][jtype])); + if (r > 0.0) fpair = -(rk/r); + else fpair = 0.0; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) evdwl = 0.5*rk*dr; + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); + } + } + } + + // bond forces + + if (atom->nmax > nmax) { + memory->sfree(s0_new); + nmax = atom->nmax; + s0_new = (double *) memory->smalloc(nmax*sizeof(double),"pair:s0_new"); + } + + // loop over my particles and their partners + // partner list contains all bond partners, so I-J appears twice + // if bond already broken, skip this partner + // first = true if this is first neighbor of particle i + + bool first; + + for (i = 0; i < nlocal; i++) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jnum = npartner[i]; + s0_new[i] = DBL_MAX; + first = true; + + for (jj = 0; jj < jnum; jj++) { + if (partner[i][jj] == 0) continue; + j = atom->map(partner[i][jj]); + + // check if lost a partner without first breaking bond + + if (j < 0) { + partner[i][jj] = 0; + continue; + } + + // compute force density, add to PD equation of motion + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + if (nonperiodic == 0) domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + delta = sqrt(cutsq[itype][jtype]); + r = sqrt(rsq); + dr = r - r0[i][jj]; + + // avoid roundoff errors + + if (fabs(dr) < 2.2204e-016) dr = 0.0; + + // scale vfrac[j] if particle j near the horizon + + if ((fabs(r0[i][jj] - delta)) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((delta - half_lc)/(2*half_lc) ) ); + else vfrac_scale = 1.0; + + stretch = dr / r0[i][jj]; + rk = (kspring[itype][jtype] * vfrac[j]) * vfrac_scale * stretch; + if (r > 0.0) fbond = -(rk/r); + else fbond = 0.0; + + f[i][0] += delx*fbond; + f[i][1] += dely*fbond; + f[i][2] += delz*fbond; + + // since I-J is double counted, set newton off & use 1/2 factor and I,I + + if (eflag) evdwl = 0.5*rk*dr; + if (evflag) ev_tally(i,i,nlocal,0, + 0.5*evdwl,0.0,0.5*fbond,delx,dely,delz); + + // find stretch in bond I-J and break if necessary + // use s0 from previous timestep + + if (stretch > MIN(s0[i],s0[j])) partner[i][jj] = 0; + + // update s0 for next timestep + + if (first) + s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch); + else + s0_new[i] = MAX(s0_new[i], + s00[itype][jtype] - (alpha[itype][jtype] * stretch)); + first = false; + } + } + + if (vflag_fdotr) virial_compute(); + + // store new s0 + + for (i = 0; i < nlocal; i++) s0[i] = s0_new[i]; +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairPeriPMB::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"); + kspring = memory->create_2d_double_array(n+1,n+1,"pair:kspring"); + s00 = memory->create_2d_double_array(n+1,n+1,"pair:s00"); + alpha = memory->create_2d_double_array(n+1,n+1,"pair:alpha"); + cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairPeriPMB::settings(int narg, char **arg) +{ + if (narg) error->all("Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairPeriPMB::coeff(int narg, char **arg) +{ + if (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 kspring_one = atof(arg[2]); + double cut_one = atof(arg[3]); + double s00_one = atof(arg[4]); + double alpha_one = atof(arg[5]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + kspring[i][j] = kspring_one; + s00[i][j] = s00_one; + alpha[i][j] = alpha_one; + cut[i][j] = cut_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 PairPeriPMB::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all("All pair coeffs are not set"); + + cutsq[i][j] = cut[i][j] * cut[i][j]; + cutsq[j][i] = cutsq[i][j]; + + // set other j,i parameters + + kspring[j][i] = kspring[i][j]; + alpha[j][i] = alpha[i][j]; + s00[j][i] = s00[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairPeriPMB::init_style() +{ + // error checks + + if (atom->map_style == 0) + error->all("Pair peri requires an atom map, see atom_modify"); + + if (atom->style_match("peri") == 0) + error->all("Pair style peri_pmb requires atom style peri"); + + if (domain->lattice == NULL) + error->all("Pair peri requires a lattice be defined"); + if (domain->lattice->xlattice != domain->lattice->ylattice || + domain->lattice->xlattice != domain->lattice->zlattice || + domain->lattice->ylattice != domain->lattice->zlattice) + error->all("Pair peri lattice is not identical in x, y, and z"); + + // if first init, create Fix needed for storing fixed neighbors + + if (ifix_peri == -1) { + char **fixarg = new char*[3]; + fixarg[0] = (char *) "PERI_NEIGH"; + fixarg[1] = (char *) "all"; + fixarg[2] = (char *) "PERI_NEIGH"; + modify->add_fix(3,fixarg); + delete [] fixarg; + } + + // find associated PERI_NEIGH fix that must exist + // could have changed locations in fix list since created + + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"PERI_NEIGH") == 0) ifix_peri = i; + + int irequest = neighbor->request(this); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairPeriPMB::write_restart(FILE *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(&kspring[i][j],sizeof(double),1,fp); + fwrite(&s00[i][j],sizeof(double),1,fp); + fwrite(&alpha[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairPeriPMB::read_restart(FILE *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(&kspring[i][j],sizeof(double),1,fp); + fread(&s00[i][j],sizeof(double),1,fp); + fread(&alpha[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&kspring[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&s00[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&alpha[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- */ + +double PairPeriPMB::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double delx0,dely0,delz0,rsq0; + double d_ij,r,dr,rk,vfrac_scale; + + double *vfrac = atom->vfrac; + double **x0 = atom->x0; + double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; + int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + + double lc = domain->lattice->xlattice; + double half_lc = 0.5*lc; + + delx0 = x0[i][0] - x0[j][0]; + dely0 = x0[i][1] - x0[j][1]; + delz0 = x0[i][2] - x0[j][2]; + if (domain->nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0); + rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; + + d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); + r = sqrt(rsq); + + double energy = 0.0; + fforce = 0.0; + + if (r < d_ij) { + dr = r - d_ij; + rk = (15.0 * kspring[itype][jtype] * vfrac[j]) * + (dr / sqrt(cutsq[itype][jtype])); + if (r > 0.0) fforce += -(rk/r); + energy += 0.5*rk*dr; + } + + int jnum = npartner[i]; + for (int jj = 0; jj < jnum; jj++) { + if (partner[i][jj] == 0) continue; + if (j < 0) continue; + if (j == atom->map(partner[i][jj])) { + dr = r - r0[i][jj]; + if (fabs(dr) < 2.2204e-016) dr = 0.0; + if ( (fabs(r0[i][jj] - sqrt(cutsq[itype][jtype]))) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((sqrt(cutsq[itype][jtype]) - half_lc)/(2*half_lc))); + else vfrac_scale = 1.0; + rk = (kspring[itype][jtype] * vfrac[j] * vfrac_scale) * + (dr / r0[i][jj]); + if (r > 0.0) fforce += -(rk/r); + energy += 0.5*rk*dr; + } + } + + return energy; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairPeriPMB::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/atom.h b/src/atom.h index f4cf022b01..8325811bff 100644 --- a/src/atom.h +++ b/src/atom.h @@ -73,6 +73,7 @@ class Atom : protected Pointers { int **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4; // per-atom array existence flags + // these can be checked before array is allocated // customize by adding new flag int molecule_flag; diff --git a/src/compute_centro_atom.h b/src/compute_centro_atom.h index 891a65a964..d09b69714b 100644 --- a/src/compute_centro_atom.h +++ b/src/compute_centro_atom.h @@ -31,8 +31,8 @@ class ComputeCentroAtom : public Compute { int nmax,maxneigh; double *distsq; int *nearest; - double *centro; class NeighList *list; + double *centro; void select(int, int, double *); void select2(int, int, double *, int *); diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp index f121f04042..657fe3f985 100644 --- a/src/compute_coord_atom.cpp +++ b/src/compute_coord_atom.cpp @@ -11,6 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include "math.h" #include "string.h" #include "stdlib.h" #include "compute_coord_atom.h" @@ -35,7 +36,8 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) : { if (narg != 4) error->all("Illegal compute coord/atom command"); - cutoff = atof(arg[3]); + double cutoff = atof(arg[3]); + cutsq = cutoff*cutoff; peratom_flag = 1; size_peratom = 0; @@ -55,7 +57,9 @@ ComputeCoordAtom::~ComputeCoordAtom() void ComputeCoordAtom::init() { - if (force->pair == NULL || cutoff > force->pair->cutforce) + if (force->pair == NULL) + error->all("Compute coord/atom requires a pair style be defined"); + if (sqrt(cutsq) > force->pair->cutforce) error->all("Compute coord/atom cutoff is longer than pairwise cutoff"); // need an occasional full neighbor list @@ -116,7 +120,6 @@ void ComputeCoordAtom::compute_peratom() double **x = atom->x; int *mask = atom->mask; int nall = atom->nlocal + atom->nghost; - double cutsq = cutoff*cutoff; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; diff --git a/src/compute_coord_atom.h b/src/compute_coord_atom.h index de1bd50a01..cd174f153a 100644 --- a/src/compute_coord_atom.h +++ b/src/compute_coord_atom.h @@ -29,9 +29,9 @@ class ComputeCoordAtom : public Compute { private: int nmax; - double cutoff; - double *coordination; + double cutsq; class NeighList *list; + double *coordination; }; } diff --git a/src/compute_erotate_sphere.cpp b/src/compute_erotate_sphere.cpp index 91fb7a3e9f..6fe43f9daa 100644 --- a/src/compute_erotate_sphere.cpp +++ b/src/compute_erotate_sphere.cpp @@ -55,8 +55,8 @@ void ComputeERotateSphere::init() if (atom->mass && !atom->shape) error->all("Compute erotate/sphere requires atom attribute shape"); - if (!atom->mass && (!atom->radius_flag || !atom->rmass_flag)) - error->all("Compute erotate/sphere requires atom attributes radius, rmass"); + if (atom->rmass && !atom->radius_flag) + error->all("Compute erotate/sphere requires atom attribute radius"); if (atom->mass) { double *mass = atom->mass; @@ -86,16 +86,16 @@ double ComputeERotateSphere::compute_scalar() double erotate = 0.0; - if (mass) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * inertia[type[i]]; - } else { + if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i]; + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * inertia[type[i]]; } MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world); diff --git a/src/compute_ke.cpp b/src/compute_ke.cpp index 7b77783c68..2db20db761 100644 --- a/src/compute_ke.cpp +++ b/src/compute_ke.cpp @@ -55,16 +55,15 @@ double ComputeKE::compute_scalar() double ke = 0.0; - if (mass) { + if (rmass) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + ke += rmass[i] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) ke += mass[type[i]] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); - - } else { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - ke += rmass[i] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); } MPI_Allreduce(&ke,&scalar,1,MPI_DOUBLE,MPI_SUM,world); diff --git a/src/compute_ke_atom.cpp b/src/compute_ke_atom.cpp index 0596e68787..432f7d5298 100644 --- a/src/compute_ke_atom.cpp +++ b/src/compute_ke_atom.cpp @@ -80,10 +80,10 @@ void ComputeKEAtom::compute_peratom() int *type = atom->type; int nlocal = atom->nlocal; - if (mass) + if (rmass) for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - ke[i] = 0.5 * mvv2e * mass[type[i]] * + ke[i] = 0.5 * mvv2e * rmass[i] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); } else ke[i] = 0.0; } @@ -91,7 +91,7 @@ void ComputeKEAtom::compute_peratom() else for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - ke[i] = 0.5 * mvv2e * rmass[i] * + ke[i] = 0.5 * mvv2e * mass[type[i]] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); } else ke[i] = 0.0; } diff --git a/src/compute_stress_atom.cpp b/src/compute_stress_atom.cpp index c10ebf453a..160f4b252e 100644 --- a/src/compute_stress_atom.cpp +++ b/src/compute_stress_atom.cpp @@ -196,10 +196,10 @@ void ComputeStressAtom::compute_peratom() int *type = atom->type; double mvv2e = force->mvv2e; - if (mass) { + if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - onemass = mvv2e * mass[type[i]]; + onemass = mvv2e * rmass[i]; stress[i][0] += onemass*v[i][0]*v[i][0]; stress[i][1] += onemass*v[i][1]*v[i][1]; stress[i][2] += onemass*v[i][2]*v[i][2]; @@ -211,7 +211,7 @@ void ComputeStressAtom::compute_peratom() } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - onemass = mvv2e * rmass[i]; + onemass = mvv2e * mass[type[i]]; stress[i][0] += onemass*v[i][0]*v[i][0]; stress[i][1] += onemass*v[i][1]*v[i][1]; stress[i][2] += onemass*v[i][2]*v[i][2]; diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp index 01c8327354..5ba16633f4 100644 --- a/src/compute_temp.cpp +++ b/src/compute_temp.cpp @@ -84,15 +84,15 @@ double ComputeTemp::compute_scalar() double t = 0.0; - if (mass) { + if (rmass) { + 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]) * rmass[i]; + } else { 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]]; - } else { - 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]) * rmass[i]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); @@ -121,8 +121,8 @@ void ComputeTemp::compute_vector() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (mass) massone = mass[type[i]]; - else massone = rmass[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp index fbbce90680..22cfab7610 100644 --- a/src/compute_temp_com.cpp +++ b/src/compute_temp_com.cpp @@ -102,12 +102,12 @@ double ComputeTempCOM::compute_scalar() vthermal[0] = v[i][0] - vbias[0]; vthermal[1] = v[i][1] - vbias[1]; vthermal[2] = v[i][2] - vbias[2]; - if (mass) - t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + - vthermal[2]*vthermal[2]) * mass[type[i]]; - else + if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; + else + t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); @@ -144,8 +144,8 @@ void ComputeTempCOM::compute_vector() vthermal[1] = v[i][1] - vbias[1]; vthermal[2] = v[i][2] - vbias[2]; - if (mass) massone = mass[type[i]]; - else massone = rmass[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp index 8fab755cce..5017dc527b 100644 --- a/src/compute_temp_deform.cpp +++ b/src/compute_temp_deform.cpp @@ -78,7 +78,8 @@ void ComputeTempDeform::init() 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 inconsistent fix deform remap option"); + error->warning("Using compute temp/deform with inconsistent " + "fix deform remap option"); break; } if (i == modify->nfix && comm->me == 0) @@ -131,12 +132,12 @@ double ComputeTempDeform::compute_scalar() vthermal[0] = v[i][0] - vstream[0]; vthermal[1] = v[i][1] - vstream[1]; vthermal[2] = v[i][2] - vstream[2]; - if (mass) - t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + - vthermal[2]*vthermal[2]) * mass[type[i]]; - else + if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; + else + t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); @@ -178,8 +179,8 @@ void ComputeTempDeform::compute_vector() vthermal[1] = v[i][1] - vstream[1]; vthermal[2] = v[i][2] - vstream[2]; - if (mass) massone = mass[type[i]]; - else massone = rmass[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp index 490a5356d1..088aca39ab 100644 --- a/src/compute_temp_partial.cpp +++ b/src/compute_temp_partial.cpp @@ -104,16 +104,16 @@ double ComputeTempPartial::compute_scalar() double t = 0.0; - if (mass) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + - zflag*v[i][2]*v[i][2]) * mass[type[i]]; - } else { + if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + zflag*v[i][2]*v[i][2]) * rmass[i]; + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + + zflag*v[i][2]*v[i][2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); @@ -142,8 +142,8 @@ void ComputeTempPartial::compute_vector() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (mass) massone = mass[type[i]]; - else massone = rmass[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; t[0] += massone * xflag*v[i][0]*v[i][0]; t[1] += massone * yflag*v[i][1]*v[i][1]; t[2] += massone * zflag*v[i][2]*v[i][2]; diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp index 395bcb81e8..5e13d3e823 100644 --- a/src/compute_temp_ramp.cpp +++ b/src/compute_temp_ramp.cpp @@ -170,12 +170,12 @@ double ComputeTempRamp::compute_scalar() vthermal[1] = v[i][1]; vthermal[2] = v[i][2]; vthermal[v_dim] -= vramp; - if (mass) - t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + - vthermal[2]*vthermal[2]) * mass[type[i]]; - else + if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; + else + t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); @@ -215,8 +215,8 @@ void ComputeTempRamp::compute_vector() vthermal[2] = v[i][2]; vthermal[v_dim] -= vramp; - if (mass) massone = mass[type[i]]; - else massone = rmass[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp index f2e06be3d6..3cdaa0606f 100644 --- a/src/compute_temp_region.cpp +++ b/src/compute_temp_region.cpp @@ -89,18 +89,18 @@ double ComputeTempRegion::compute_scalar() int count = 0; double t = 0.0; - if (mass) { + if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { count++; - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - mass[type[i]]; + t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { count++; - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; + t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + mass[type[i]]; } } @@ -136,8 +136,8 @@ void ComputeTempRegion::compute_vector() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { - if (mass) massone = mass[type[i]]; - else massone = rmass[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp index d338e78fb6..31bff920e9 100644 --- a/src/compute_temp_sphere.cpp +++ b/src/compute_temp_sphere.cpp @@ -152,7 +152,14 @@ double ComputeTempSphere::compute_scalar() double t = 0.0; - if (mass) { + if (rmass) { + 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]) * rmass[i]; + t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * INERTIA*radius[i]*radius[i]*rmass[i]; + } + } else { 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]) * @@ -160,13 +167,6 @@ double ComputeTempSphere::compute_scalar() t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + omega[i][2]*omega[i][2]) * inertia[type[i]]; } - } else { - 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]) * rmass[i]; - t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * INERTIA*radius[i]*radius[i]*rmass[i]; - } } if (tempbias) tbias->restore_bias_all(); @@ -202,26 +202,7 @@ void ComputeTempSphere::compute_vector() double massone,inertiaone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; - if (mass) { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - massone = mass[type[i]]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; - - inertiaone = inertia[type[i]]; - t[0] += inertiaone * omega[i][0]*omega[i][0]; - t[1] += inertiaone * omega[i][1]*omega[i][1]; - t[2] += inertiaone * omega[i][2]*omega[i][2]; - t[3] += inertiaone * omega[i][0]*omega[i][1]; - t[4] += inertiaone * omega[i][0]*omega[i][2]; - t[5] += inertiaone * omega[i][1]*omega[i][2]; - } - } else { + if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { massone = rmass[i]; @@ -240,6 +221,25 @@ void ComputeTempSphere::compute_vector() t[4] += inertiaone * omega[i][0]*omega[i][2]; t[5] += inertiaone * omega[i][1]*omega[i][2]; } + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + massone = mass[type[i]]; + t[0] += massone * v[i][0]*v[i][0]; + t[1] += massone * v[i][1]*v[i][1]; + t[2] += massone * v[i][2]*v[i][2]; + t[3] += massone * v[i][0]*v[i][1]; + t[4] += massone * v[i][0]*v[i][2]; + t[5] += massone * v[i][1]*v[i][2]; + + inertiaone = inertia[type[i]]; + t[0] += inertiaone * omega[i][0]*omega[i][0]; + t[1] += inertiaone * omega[i][1]*omega[i][1]; + t[2] += inertiaone * omega[i][2]*omega[i][2]; + t[3] += inertiaone * omega[i][0]*omega[i][1]; + t[4] += inertiaone * omega[i][0]*omega[i][2]; + t[5] += inertiaone * omega[i][1]*omega[i][2]; + } } if (tempbias) tbias->restore_bias_all(); diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 09f6753c7f..ea6fecb596 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -343,15 +343,15 @@ int DumpCustom::count() ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == MASS) { - if (atom->mass) { + if (atom->rmass) { + ptr = atom->rmass; + nstride = 1; + } else { double *mass = atom->mass; int *type = atom->type; for (i = 0; i < nlocal; i++) dchoose[i] = mass[type[i]]; ptr = dchoose; nstride = 1; - } else { - ptr = atom->rmass; - nstride = 1; } } else if (thresh_array[ithresh] == X) { @@ -1450,16 +1450,16 @@ void DumpCustom::pack_mass(int n) double *rmass = atom->rmass; int nlocal = atom->nlocal; - if (mass) { + if (rmass) { for (int i = 0; i < nlocal; i++) if (choose[i]) { - buf[n] = mass[type[i]]; + buf[n] = rmass[i]; n += size_one; } } else { for (int i = 0; i < nlocal; i++) if (choose[i]) { - buf[n] = rmass[i]; + buf[n] = mass[type[i]]; n += size_one; } } diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp index fd9c36cdad..dfdafd96f1 100644 --- a/src/fix_ave_spatial.cpp +++ b/src/fix_ave_spatial.cpp @@ -615,8 +615,8 @@ void FixAveSpatial::end_of_step() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) - if (mass) values_one[layer[i]][m] += mass[type[i]]; - else values_one[layer[i]][m] += rmass[i]; + if (rmass) values_one[layer[i]][m] += rmass[i]; + else values_one[layer[i]][m] += mass[type[i]]; // COMPUTE adds its scalar or vector component to values // invoke compute if not previously invoked diff --git a/src/fix_dt_reset.cpp b/src/fix_dt_reset.cpp index 0fe6f0f101..b1275c65aa 100644 --- a/src/fix_dt_reset.cpp +++ b/src/fix_dt_reset.cpp @@ -152,8 +152,8 @@ void FixDtReset::end_of_step() if (mask[i] & groupbit) { vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; fsq = f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2]; - if (mass) ms = mass[type[i]]; - else ms = rmass[i]; + if (rmass) ms = rmass[i]; + else ms = mass[type[i]]; asq = fsq/ms/ms; bound[0] = MAX(bound[0],vsq); bound[1] = MAX(bound[1],asq); diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp index 9dda46ea30..17b0f5bd51 100644 --- a/src/fix_gravity.cpp +++ b/src/fix_gravity.cpp @@ -152,10 +152,10 @@ void FixGravity::post_force(int vflag) int nlocal = atom->nlocal; double massone; - if (mass) { + if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - massone = mass[type[i]]; + massone = rmass[i]; f[i][0] += massone*xacc; f[i][1] += massone*yacc; f[i][2] += massone*zacc; @@ -163,9 +163,10 @@ void FixGravity::post_force(int vflag) } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - f[i][0] += rmass[i]*xacc; - f[i][1] += rmass[i]*yacc; - f[i][2] += rmass[i]*zacc; + massone = mass[type[i]]; + f[i][0] += massone*xacc; + f[i][1] += massone*yacc; + f[i][2] += massone*zacc; } } } diff --git a/src/fix_line_force.cpp b/src/fix_line_force.cpp index 530556b380..389b80f2dc 100644 --- a/src/fix_line_force.cpp +++ b/src/fix_line_force.cpp @@ -105,5 +105,3 @@ void FixLineForce::min_post_force(int vflag) { post_force(vflag); } - - diff --git a/src/fix_nve.cpp b/src/fix_nve.cpp index ec1258b5f9..f8789c276e 100644 --- a/src/fix_nve.cpp +++ b/src/fix_nve.cpp @@ -74,10 +74,10 @@ void FixNVE::initial_integrate(int vflag) int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (mass) { + if (rmass) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; + dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; @@ -90,7 +90,7 @@ void FixNVE::initial_integrate(int vflag) } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; + 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]; @@ -117,10 +117,10 @@ void FixNVE::final_integrate() int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (mass) { + if (rmass) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; + dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; @@ -130,7 +130,7 @@ void FixNVE::final_integrate() } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; + 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]; diff --git a/src/fix_nve_limit.cpp b/src/fix_nve_limit.cpp index 3cbd46eb30..c4465487cd 100644 --- a/src/fix_nve_limit.cpp +++ b/src/fix_nve_limit.cpp @@ -84,10 +84,10 @@ void FixNVELimit::initial_integrate(int vflag) int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (mass) { + if (rmass) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; + dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; @@ -110,7 +110,7 @@ void FixNVELimit::initial_integrate(int vflag) } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; + 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]; @@ -147,10 +147,10 @@ void FixNVELimit::final_integrate() int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (mass) { + if (rmass) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; + dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; @@ -169,7 +169,7 @@ void FixNVELimit::final_integrate() } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; + 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]; diff --git a/src/fix_nve_sphere.cpp b/src/fix_nve_sphere.cpp index 41650c3ad3..aced7ec097 100644 --- a/src/fix_nve_sphere.cpp +++ b/src/fix_nve_sphere.cpp @@ -92,8 +92,8 @@ void FixNVESphere::init() if (atom->mass && !atom->shape) error->all("Fix nve/sphere requires atom attribute shape"); - if (!atom->mass && (!atom->radius_flag || !atom->rmass_flag)) - error->all("Fix nve/sphere requires atom attributes radius, rmass"); + if (atom->rmass && !atom->radius_flag) + error->all("Fix nve/sphere requires atom attribute radius"); if (atom->mass) { double **shape = atom->shape; @@ -137,7 +137,25 @@ void FixNVESphere::initial_integrate(int vflag) // update v,x,omega for all particles // d_omega/dt = torque / inertia - if (mass) { + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[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]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } + } + + } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { itype = type[i]; @@ -155,24 +173,6 @@ void FixNVESphere::initial_integrate(int vflag) omega[i][2] += dtirotate * torque[i][2]; } } - - } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; - } - } } // update mu for dipoles @@ -228,7 +228,23 @@ void FixNVESphere::final_integrate() dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]); } - if (mass) { + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } + } + + } else { + dtfrotate = dtf / INERTIA; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { itype = type[i]; @@ -243,21 +259,5 @@ void FixNVESphere::final_integrate() omega[i][2] += dtirotate * torque[i][2]; } } - - } else { - dtfrotate = dtf / INERTIA; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; - } - } } } diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp index a55c828daf..e13e04c539 100644 --- a/src/fix_shake.cpp +++ b/src/fix_shake.cpp @@ -121,10 +121,10 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : type_flag[i] = 1; } else if (mode == 'm') { - double rmass = atof(arg[next]); - if (rmass == 0.0) error->all("Invalid atom mass for fix shake"); + double massone = atof(arg[next]); + if (massone == 0.0) error->all("Invalid atom mass for fix shake"); if (nmass == atom->ntypes) error->all("Too many masses for fix shake"); - mass_list[nmass++] = rmass; + mass_list[nmass++] = massone; } else error->all("Illegal fix shake command"); next++; @@ -428,6 +428,7 @@ void FixShake::pre_neighbor() v = atom->v; f = atom->f; mass = atom->mass; + rmass = atom->rmass; type = atom->type; nlocal = atom->nlocal; @@ -562,7 +563,7 @@ void FixShake::find_clusters() { int i,j,m,n; int flag,flag_all,messtag,loop,nbuf,nbufmax,size; - double imass,jmass,rmass; + double massone; int *buf,*bufcopy; MPI_Request request; MPI_Status status; @@ -575,6 +576,7 @@ void FixShake::find_clusters() int *type = atom->type; int *mask = atom->mask; double *mass = atom->mass; + double *rmass = atom->rmass; int **bond_type = atom->bond_type; int **angle_type = atom->angle_type; int **nspecial = atom->nspecial; @@ -596,6 +598,7 @@ void FixShake::find_clusters() // partner_tag[i][] = global IDs of each partner // partner_mask[i][] = mask of each partner // partner_type[i][] = type of each partner + // partner_massflag[i][] = 1 if partner meets mass criterion, 0 if not // partner_bondtype[i][] = type of bond attached to each partner // partner_shake[i][] = 1 if SHAKE bonded to partner, 0 if not // partner_nshake[i][] = nshake value for each partner @@ -614,6 +617,8 @@ void FixShake::find_clusters() memory->create_2d_int_array(nlocal,max,"shake:partner_mask"); int **partner_type = memory->create_2d_int_array(nlocal,max,"shake:partner_type"); + int **partner_massflag = + memory->create_2d_int_array(nlocal,max,"shake:partner_massflag"); int **partner_bondtype = memory->create_2d_int_array(nlocal,max,"shake:partner_bondtype"); int **partner_shake = @@ -631,33 +636,42 @@ void FixShake::find_clusters() } // ----------------------------------------------------- - // set partner_mask, partner_type, partner_bondtype for bonded partners + // set partner_mask, partner_type, partner_massflag, partner_bondtype + // for bonded partners // requires communication for off-proc partners // ----------------------------------------------------- - // fill in mask, type, bondtype if own bond partner - // info to store in buf for each off-proc bond = - // 2 atoms IDs in bond, space for mask, type, bondtype + // fill in mask, type, massflag, bondtype if own bond partner + // info to store in buf for each off-proc bond = nper = 6 + // 2 atoms IDs in bond, space for mask, type, massflag, bondtype // nbufmax = largest buffer needed to hold info from any proc + int nper = 6; + nbuf = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < npartner[i]; j++) { partner_mask[i][j] = 0; partner_type[i][j] = 0; + partner_massflag[i][j] = 0; partner_bondtype[i][j] = 0; m = atom->map(partner_tag[i][j]); if (m >= 0 && m < nlocal) { partner_mask[i][j] = mask[m]; partner_type[i][j] = type[m]; + if (nmass) { + if (rmass) massone = rmass[m]; + else massone = mass[type[m]]; + partner_massflag[i][j] = masscheck(massone); + } n = bondfind(i,tag[i],partner_tag[i][j]); if (n >= 0) partner_bondtype[i][j] = bond_type[i][n]; else { n = bondfind(m,tag[i],partner_tag[i][j]); if (n >= 0) partner_bondtype[i][j] = bond_type[m][n]; } - } else nbuf += 5; + } else nbuf += nper; } } @@ -677,18 +691,20 @@ void FixShake::find_clusters() buf[size+1] = partner_tag[i][j]; buf[size+2] = 0; buf[size+3] = 0; + buf[size+4] = 0; n = bondfind(i,tag[i],partner_tag[i][j]); - if (n >= 0) buf[size+4] = bond_type[i][n]; - else buf[size+4] = 0; - size += 5; + if (n >= 0) buf[size+5] = bond_type[i][n]; + else buf[size+5] = 0; + size += nper; } } } // cycle buffer around ring of procs back to self // when receive buffer, scan bond partner IDs for atoms I own - // if I own partner, fill in mask and type, search for bond with 1st atom - // and fill in bondtype + // if I own partner: + // fill in mask and type and massflag + // search for bond with 1st atom and fill in bondtype messtag = 1; for (loop = 0; loop < nprocs; loop++) { @@ -698,12 +714,17 @@ void FixShake::find_clusters() if (m >= 0 && m < nlocal) { buf[i+2] = mask[m]; buf[i+3] = type[m]; - if (buf[i+4] == 0) { + if (nmass) { + if (rmass) massone = rmass[m]; + else massone = mass[type[m]]; + buf[i+4] = masscheck(massone); + } + if (buf[i+5] == 0) { n = bondfind(m,buf[i],buf[i+1]); - if (n >= 0) buf[i+4] = bond_type[m][n]; + if (n >= 0) buf[i+5] = bond_type[m][n]; } } - i += 5; + i += nper; } if (me != next) { MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); @@ -723,8 +744,9 @@ void FixShake::find_clusters() if (buf[m+1] == partner_tag[i][j]) break; partner_mask[i][j] = buf[m+2]; partner_type[i][j] = buf[m+3]; - partner_bondtype[i][j] = buf[m+4]; - m += 5; + partner_massflag[i][j] = buf[m+4]; + partner_bondtype[i][j] = buf[m+5]; + m += nper; } delete [] buf; @@ -780,15 +802,17 @@ void FixShake::find_clusters() continue; } if (nmass) { - imass = mass[type[i]]; - jmass = mass[partner_type[i][j]]; - for (m = 0; m < nmass; m++) { - rmass = mass_list[m]; - if (fabs(rmass-imass) <= MASSDELTA || - fabs(rmass-jmass) <= MASSDELTA) { + if (partner_massflag[i][j]) { + partner_shake[i][j] = 1; + nshake[i]++; + continue; + } else { + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + if (masscheck(massone)) { partner_shake[i][j] = 1; nshake[i]++; - break; + continue; } } } @@ -1049,6 +1073,7 @@ void FixShake::find_clusters() memory->destroy_2d_int_array(partner_tag); memory->destroy_2d_int_array(partner_mask); memory->destroy_2d_int_array(partner_type); + memory->destroy_2d_int_array(partner_massflag); memory->destroy_2d_int_array(partner_bondtype); memory->destroy_2d_int_array(partner_shake); memory->destroy_2d_int_array(partner_nshake); @@ -1124,6 +1149,18 @@ void FixShake::find_clusters() } } +/* ---------------------------------------------------------------------- + check if massone is within MASSDELTA of any mass in mass_list + return 1 if yes, 0 if not +------------------------------------------------------------------------- */ + +int FixShake::masscheck(double massone) +{ + for (int i = 0; i < nmass; i++) + if (fabs(mass_list[i]-massone) <= MASSDELTA) return 1; + return 0; +} + /* ---------------------------------------------------------------------- update the unconstrained position of each atom only for SHAKE clusters, else set to 0.0 @@ -1134,13 +1171,24 @@ void FixShake::unconstrained_update() { double dtfmsq; - for (int i = 0; i < nlocal; i++) { - if (shake_flag[i]) { - dtfmsq = dtfsq / mass[type[i]]; - xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; - xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; - xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; - } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (shake_flag[i]) { + dtfmsq = dtfsq / rmass[i]; + xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; + xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; + xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; + } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; + } + } else { + for (int i = 0; i < nlocal; i++) { + if (shake_flag[i]) { + dtfmsq = dtfsq / mass[type[i]]; + xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; + xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; + xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; + } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; + } } } @@ -1150,6 +1198,7 @@ void FixShake::shake2(int m) { int nlist,list[2]; double v[6]; + double invmass0,invmass1; // local atom IDs and constraint distances @@ -1180,8 +1229,13 @@ void FixShake::shake2(int m) // a,b,c = coeffs in quadratic equation for lamda - double invmass0 = 1.0/mass[type[i0]]; - double invmass1 = 1.0/mass[type[i1]]; + if (rmass) { + invmass0 = 1.0/rmass[i0]; + invmass1 = 1.0/rmass[i1]; + } else { + invmass0 = 1.0/mass[type[i0]]; + invmass1 = 1.0/mass[type[i1]]; + } double a = (invmass0+invmass1)*(invmass0+invmass1) * r01sq; double b = 2.0 * (invmass0+invmass1) * @@ -1243,6 +1297,7 @@ void FixShake::shake3(int m) { int nlist,list[3]; double v[6]; + double invmass0,invmass1,invmass2; // local atom IDs and constraint distances @@ -1289,9 +1344,15 @@ void FixShake::shake3(int m) // matrix coeffs and rhs for lamda equations - double invmass0 = 1.0/mass[type[i0]]; - double invmass1 = 1.0/mass[type[i1]]; - double invmass2 = 1.0/mass[type[i2]]; + if (rmass) { + invmass0 = 1.0/rmass[i0]; + invmass1 = 1.0/rmass[i1]; + invmass2 = 1.0/rmass[i2]; + } else { + invmass0 = 1.0/mass[type[i0]]; + invmass1 = 1.0/mass[type[i1]]; + invmass2 = 1.0/mass[type[i2]]; + } double a11 = 2.0 * (invmass0+invmass1) * (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); @@ -1401,6 +1462,7 @@ void FixShake::shake4(int m) { int nlist,list[4]; double v[6]; + double invmass0,invmass1,invmass2,invmass3; // local atom IDs and constraint distances @@ -1463,10 +1525,17 @@ void FixShake::shake4(int m) // matrix coeffs and rhs for lamda equations - double invmass0 = 1.0/mass[type[i0]]; - double invmass1 = 1.0/mass[type[i1]]; - double invmass2 = 1.0/mass[type[i2]]; - double invmass3 = 1.0/mass[type[i3]]; + if (rmass) { + invmass0 = 1.0/rmass[i0]; + invmass1 = 1.0/rmass[i1]; + invmass2 = 1.0/rmass[i2]; + invmass3 = 1.0/rmass[i3]; + } else { + invmass0 = 1.0/mass[type[i0]]; + invmass1 = 1.0/mass[type[i1]]; + invmass2 = 1.0/mass[type[i2]]; + invmass3 = 1.0/mass[type[i3]]; + } double a11 = 2.0 * (invmass0+invmass1) * (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); @@ -1636,6 +1705,7 @@ void FixShake::shake3angle(int m) { int nlist,list[3]; double v[6]; + double invmass0,invmass1,invmass2; // local atom IDs and constraint distances @@ -1697,9 +1767,15 @@ void FixShake::shake3angle(int m) // matrix coeffs and rhs for lamda equations - double invmass0 = 1.0/mass[type[i0]]; - double invmass1 = 1.0/mass[type[i1]]; - double invmass2 = 1.0/mass[type[i2]]; + if (rmass) { + invmass0 = 1.0/rmass[i0]; + invmass1 = 1.0/rmass[i1]; + invmass2 = 1.0/rmass[i2]; + } else { + invmass0 = 1.0/mass[type[i0]]; + invmass1 = 1.0/mass[type[i1]]; + invmass2 = 1.0/mass[type[i2]]; + } double a11 = 2.0 * (invmass0+invmass1) * (s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]); @@ -2204,20 +2280,39 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop) double invmass,dtfmsq; int jlevel; - for (int i = 0; i < nlocal; i++) { - if (shake_flag[i]) { - invmass = 1.0 / mass[type[i]]; - dtfmsq = dtfsq * invmass; - xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; - xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; - xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; - for (jlevel = 0; jlevel < ilevel; jlevel++) { - dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; - xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; - xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; - xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; - } - } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (shake_flag[i]) { + invmass = 1.0 / rmass[i]; + dtfmsq = dtfsq * invmass; + xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; + xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; + xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; + for (jlevel = 0; jlevel < ilevel; jlevel++) { + dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; + xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; + xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; + xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; + } + } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; + } + + } else { + for (int i = 0; i < nlocal; i++) { + if (shake_flag[i]) { + invmass = 1.0 / mass[type[i]]; + dtfmsq = dtfsq * invmass; + xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; + xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; + xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; + for (jlevel = 0; jlevel < ilevel; jlevel++) { + dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; + xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; + xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; + xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; + } + } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; + } } // communicate results if necessary diff --git a/src/fix_shake.h b/src/fix_shake.h index 3512bc43ea..74aa51ebbb 100644 --- a/src/fix_shake.h +++ b/src/fix_shake.h @@ -52,7 +52,7 @@ class FixShake : public Fix { int *bond_flag,*angle_flag; // bond/angle types to constrain int *type_flag; // constrain bonds to these types double *mass_list; // constrain bonds to these masses - int nmass; + int nmass; // # of masses in mass_list double *bond_distance,*angle_distance; // constraint distances @@ -62,7 +62,7 @@ class FixShake : public Fix { double *step_respa; double **x,**v,**f; // local ptrs to atom class quantities - double *mass; + double *mass,*rmass; int *type; int nlocal; // atom-based arrays @@ -93,6 +93,7 @@ class FixShake : public Fix { double *a_ave_all,*a_max_all,*a_min_all; void find_clusters(); + int masscheck(double); void unconstrained_update(); void shake2(int); void shake3(int); diff --git a/src/fix_thermal_conductivity.cpp b/src/fix_thermal_conductivity.cpp index 6bd2523fb2..7321d66402 100644 --- a/src/fix_thermal_conductivity.cpp +++ b/src/fix_thermal_conductivity.cpp @@ -182,8 +182,8 @@ void FixThermalConductivity::end_of_step() if (coord >= slablo_lo && coord < slablo_hi) { ke = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; - if (mass) ke *= 0.5*mass[type[i]]; - else ke *= 0.5*rmass[i]; + if (rmass) ke *= 0.5*rmass[i]; + else ke *= 0.5*mass[type[i]]; if (nlo < nswap || ke > ke_lo[nswap-1]) { for (insert = nlo-1; insert >= 0; insert--) if (ke < ke_lo[insert]) break; @@ -200,8 +200,8 @@ void FixThermalConductivity::end_of_step() if (coord >= slabhi_lo && coord < slabhi_hi) { ke = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; - if (mass) ke *= 0.5*mass[type[i]]; - else ke *= 0.5*rmass[i]; + if (rmass) ke *= 0.5*rmass[i]; + else ke *= 0.5*mass[type[i]]; if (nhi < nswap || ke < ke_hi[nswap-1]) { for (insert = nhi-1; insert >= 0; insert--) if (ke > ke_hi[insert]) break; diff --git a/src/fix_viscosity.cpp b/src/fix_viscosity.cpp index 81085eeb2c..d6dcc1fc77 100644 --- a/src/fix_viscosity.cpp +++ b/src/fix_viscosity.cpp @@ -176,8 +176,6 @@ void FixViscosity::end_of_step() 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; @@ -232,6 +230,9 @@ void FixViscosity::end_of_step() // exchange momenta between the 2 particles // if I own both particles just swap, else point2point comm of vel,mass + double *mass = atom->mass; + double *rmass = atom->rmass; + int ipos,ineg; double sbuf[2],rbuf[2]; @@ -254,10 +255,10 @@ void FixViscosity::end_of_step() ipos = pos_index[ipositive++]; ineg = neg_index[inegative++]; rbuf[0] = v[ipos][vdim]; - if (mass) rbuf[1] = mass[type[ipos]]; - else rbuf[1] = rmass[ipos]; + if (rmass) rbuf[1] = rmass[ipos]; + else rbuf[1] = mass[type[ipos]]; sbuf[0] = v[ineg][vdim]; - if (mass) sbuf[1] = mass[type[ineg]]; + if (rmass) sbuf[1] = mass[type[ineg]]; else sbuf[1] = rmass[ineg]; v[ineg][vdim] = rbuf[0] * rbuf[1]/sbuf[1]; v[ipos][vdim] = sbuf[0] * sbuf[1]/rbuf[1]; @@ -266,8 +267,8 @@ void FixViscosity::end_of_step() } else if (me == all[0].proc) { ipos = pos_index[ipositive++]; sbuf[0] = v[ipos][vdim]; - if (mass) sbuf[1] = mass[type[ipos]]; - else sbuf[1] = rmass[ipos]; + if (rmass) sbuf[1] = rmass[ipos]; + else sbuf[1] = mass[type[ipos]]; MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[1].proc,0, rbuf,2,MPI_DOUBLE,all[1].proc,0,world,&status); v[ipos][vdim] = rbuf[0] * rbuf[1]/sbuf[1]; @@ -276,8 +277,8 @@ void FixViscosity::end_of_step() } else if (me == all[1].proc) { ineg = neg_index[inegative++]; sbuf[0] = v[ineg][vdim]; - if (mass) sbuf[1] = mass[type[ineg]]; - else sbuf[1] = rmass[ineg]; + if (rmass) sbuf[1] = rmass[ineg]; + else sbuf[1] = mass[type[ineg]]; MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[0].proc,0, rbuf,2,MPI_DOUBLE,all[0].proc,0,world,&status); v[ineg][vdim] = rbuf[0] * rbuf[1]/sbuf[1]; diff --git a/src/group.cpp b/src/group.cpp index 39b881cbdf..42b392be4b 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -446,12 +446,12 @@ double Group::mass(int igroup) double one = 0.0; - if (mass) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) one += mass[type[i]]; - } else { + if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) one += rmass[i]; + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) one += mass[type[i]]; } double all; @@ -550,13 +550,13 @@ void Group::xcm(int igroup, double masstotal, double *cm) int xbox,ybox,zbox; double massone; - if (mass) { + if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { xbox = (image[i] & 1023) - 512; ybox = (image[i] >> 10 & 1023) - 512; zbox = (image[i] >> 20) - 512; - massone = mass[type[i]]; + massone = rmass[i]; cmone[0] += (x[i][0] + xbox*xprd) * massone; cmone[1] += (x[i][1] + ybox*yprd) * massone; cmone[2] += (x[i][2] + zbox*zprd) * massone; @@ -567,7 +567,7 @@ void Group::xcm(int igroup, double masstotal, double *cm) xbox = (image[i] & 1023) - 512; ybox = (image[i] >> 10 & 1023) - 512; zbox = (image[i] >> 20) - 512; - massone = rmass[i]; + massone = mass[type[i]]; cmone[0] += (x[i][0] + xbox*xprd) * massone; cmone[1] += (x[i][1] + ybox*yprd) * massone; cmone[2] += (x[i][2] + zbox*zprd) * massone; @@ -599,10 +599,10 @@ void Group::vcm(int igroup, double masstotal, double *cm) double p[3],massone; p[0] = p[1] = p[2] = 0.0; - if (mass) { + if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - massone = mass[type[i]]; + massone = rmass[i]; p[0] += v[i][0]*massone; p[1] += v[i][1]*massone; p[2] += v[i][2]*massone; @@ -610,7 +610,7 @@ void Group::vcm(int igroup, double masstotal, double *cm) } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - massone = rmass[i]; + massone = mass[type[i]]; p[0] += v[i][0]*massone; p[1] += v[i][1]*massone; p[2] += v[i][2]*massone; @@ -665,16 +665,16 @@ double Group::ke(int igroup) double one = 0.0; - if (mass) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - mass[type[i]]; - } else { + if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + mass[type[i]]; } double all; @@ -715,8 +715,8 @@ double Group::gyration(int igroup, double masstotal, double *cm) dx = (x[i][0] + xbox*xprd) - cm[0]; dy = (x[i][1] + ybox*yprd) - cm[1]; dz = (x[i][2] + zbox*zprd) - cm[2]; - if (mass) massone = mass[type[i]]; - else massone = rmass[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; rg += (dx*dx + dy*dy + dz*dz) * massone; } double rg_all; @@ -759,8 +759,8 @@ void Group::angmom(int igroup, double *cm, double *lmom) dx = (x[i][0] + xbox*xprd) - cm[0]; dy = (x[i][1] + ybox*yprd) - cm[1]; dz = (x[i][2] + zbox*zprd) - cm[2]; - if (mass) massone = mass[type[i]]; - else massone = rmass[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; p[0] += massone * (dy*v[i][2] - dz*v[i][1]); p[1] += massone * (dz*v[i][0] - dx*v[i][2]); p[2] += massone * (dx*v[i][1] - dy*v[i][0]); @@ -806,8 +806,8 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3]) dx = (x[i][0] + xbox*xprd) - cm[0]; dy = (x[i][1] + ybox*yprd) - cm[1]; dz = (x[i][2] + zbox*zprd) - cm[2]; - if (mass) massone = mass[type[i]]; - else massone = rmass[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; ione[0][0] += massone * (dy*dy + dz*dz); ione[1][1] += massone * (dx*dx + dz*dz); ione[2][2] += massone * (dx*dx + dy*dy); diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index d1180ca501..f3fdffbe14 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -714,11 +714,3 @@ double PairLJCut::single(int i, int j, int itype, int jtype, double rsq, offset[itype][jtype]; return factor_lj*philj; } - -/* ---------------------------------------------------------------------- */ - -void *PairLJCut::extract(char *str) -{ - if (strcmp(str,"diameter") == 0) return (void *) sigma; - return NULL; -} diff --git a/src/pair_lj_cut.h b/src/pair_lj_cut.h index 715e9f194c..ab6f7d0ceb 100644 --- a/src/pair_lj_cut.h +++ b/src/pair_lj_cut.h @@ -37,7 +37,6 @@ class PairLJCut : public Pair { void compute_inner(); void compute_middle(); void compute_outer(int, int); - void *extract(char *); protected: double cut_global; diff --git a/src/velocity.cpp b/src/velocity.cpp index 7ba5361f39..70ca09065b 100644 --- a/src/velocity.cpp +++ b/src/velocity.cpp @@ -225,8 +225,8 @@ void Velocity::create(int narg, char **arg) m = atom->map(i); if (m >= 0 && m < nlocal) { if (mask[m] & groupbit) { - if (mass) factor = 1.0/sqrt(mass[type[m]]); - else factor = 1.0/sqrt(rmass[m]); + if (rmass) factor = 1.0/sqrt(rmass[m]); + factor = 1.0/sqrt(mass[type[m]]); v[m][0] = vx * factor; v[m][1] = vy * factor; if (dimension == 3) v[m][2] = vz * factor; @@ -257,8 +257,8 @@ void Velocity::create(int narg, char **arg) vy = random->gaussian(); vz = random->gaussian(); } - if (mass) factor = 1.0/sqrt(mass[type[i]]); - else factor = 1.0/sqrt(rmass[i]); + if (rmass) factor = 1.0/sqrt(rmass[i]); + factor = 1.0/sqrt(mass[type[i]]); v[i][0] = vx * factor; v[i][1] = vy * factor; if (dimension == 3) v[i][2] = vz * factor; @@ -283,8 +283,8 @@ void Velocity::create(int narg, char **arg) vz = random->gaussian(); } - if (mass) factor = 1.0/sqrt(mass[type[i]]); - else factor = 1.0/sqrt(rmass[i]); + if (rmass) factor = 1.0/sqrt(rmass[i]); + factor = 1.0/sqrt(mass[type[i]]); v[i][0] = vx * factor; v[i][1] = vy * factor; if (dimension == 3) v[i][2] = vz * factor;