diff --git a/src/DIPOLE/atom_vec_dipole.cpp b/src/DIPOLE/atom_vec_dipole.cpp index 2b94516b7e..25377ae828 100644 --- a/src/DIPOLE/atom_vec_dipole.cpp +++ b/src/DIPOLE/atom_vec_dipole.cpp @@ -737,7 +737,7 @@ void AtomVecDipole::create_atom(int itype, double *coord) x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; mask[nlocal] = 1; - image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | + image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMAX; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; diff --git a/src/GPU/pair_eam_gpu.cpp b/src/GPU/pair_eam_gpu.cpp index 113a9df059..25e6ba8140 100644 --- a/src/GPU/pair_eam_gpu.cpp +++ b/src/GPU/pair_eam_gpu.cpp @@ -5,7 +5,7 @@ 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 + 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. @@ -38,30 +38,30 @@ using namespace LAMMPS_NS; // External functions from cuda library for atom decomposition int eam_gpu_init(const int ntypes, double host_cutforcesq, - int **host_type2rhor, int **host_type2z2r, + int **host_type2rhor, int **host_type2z2r, int *host_type2frho, double ***host_rhor_spline, - double ***host_z2r_spline, double ***host_frho_spline, + double ***host_z2r_spline, double ***host_frho_spline, double rdr, double rdrho, int nrhor, int nrho, int nz2r, - int nfrho, int nr, const int nlocal, const int nall, - const int max_nbors, const int maxspecial, - const double cell_size, int &gpu_mode, FILE *screen, - int &fp_size); + int nfrho, int nr, const int nlocal, const int nall, + const int max_nbors, const int maxspecial, + const double cell_size, int &gpu_mode, FILE *screen, + int &fp_size); void eam_gpu_clear(); int** eam_gpu_compute_n(const int ago, const int inum_full, const int nall, - double **host_x, int *host_type, double *sublo, - double *subhi, int *tag, int **nspecial, int **special, - const bool eflag, const bool vflag, const bool eatom, - const bool vatom, int &host_start, int **ilist, - int **jnum, const double cpu_time, bool &success, - int &inum, void **fp_ptr); -void eam_gpu_compute(const int ago, const int inum_full, const int nlocal, - const int nall,double **host_x, int *host_type, - int *ilist, int *numj, int **firstneigh, - const bool eflag, const bool vflag, - const bool eatom, const bool vatom, int &host_start, - const double cpu_time, bool &success, void **fp_ptr); + double **host_x, int *host_type, double *sublo, + double *subhi, int *tag, int **nspecial, int **special, + const bool eflag, const bool vflag, const bool eatom, + const bool vatom, int &host_start, int **ilist, + int **jnum, const double cpu_time, bool &success, + int &inum, void **fp_ptr); +void eam_gpu_compute(const int ago, const int inum_full, const int nlocal, + const int nall,double **host_x, int *host_type, + int *ilist, int *numj, int **firstneigh, + const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + const double cpu_time, bool &success, void **fp_ptr); void eam_gpu_compute_force(int *ilist, const bool eflag, const bool vflag, - const bool eatom, const bool vatom); + const bool eatom, const bool vatom); double eam_gpu_bytes(); /* ---------------------------------------------------------------------- */ @@ -70,7 +70,7 @@ PairEAMGPU::PairEAMGPU(LAMMPS *lmp) : PairEAM(lmp), gpu_mode(GPU_FORCE) { respa_enable = 0; cpu_time = 0.0; - GPU_EXTRA::gpu_ready(lmp->modify, lmp->error); + GPU_EXTRA::gpu_ready(lmp->modify, lmp->error); } /* ---------------------------------------------------------------------- @@ -100,44 +100,44 @@ void PairEAMGPU::compute(int eflag, int vflag) evdwl = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = eflag_global = eflag_atom = 0; - + int nlocal = atom->nlocal; int newton_pair = force->newton_pair; // compute density on each atom on GPU - int nall = atom->nlocal + atom->nghost; + int nall = atom->nlocal + atom->nghost; int inum, host_start, inum_dev; - + bool success = true; - int *ilist, *numneigh, **firstneigh; - if (gpu_mode != GPU_FORCE) { + int *ilist, *numneigh, **firstneigh; + if (gpu_mode != GPU_FORCE) { inum = atom->nlocal; firstneigh = eam_gpu_compute_n(neighbor->ago, inum, nall, atom->x, - atom->type, domain->sublo, domain->subhi, - atom->tag, atom->nspecial, atom->special, - eflag, vflag, eflag_atom, vflag_atom, - host_start, &ilist, &numneigh, cpu_time, - success, inum_dev, &fp_pinned); + atom->type, domain->sublo, domain->subhi, + atom->tag, atom->nspecial, atom->special, + eflag, vflag, eflag_atom, vflag_atom, + host_start, &ilist, &numneigh, cpu_time, + success, inum_dev, &fp_pinned); } else { // gpu_mode == GPU_FORCE inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; eam_gpu_compute(neighbor->ago, inum, nlocal, nall, atom->x, atom->type, - ilist, numneigh, firstneigh, eflag, vflag, eflag_atom, - vflag_atom, host_start, cpu_time, success, &fp_pinned); + ilist, numneigh, firstneigh, eflag, vflag, eflag_atom, + vflag_atom, host_start, cpu_time, success, &fp_pinned); } - + if (!success) error->one(FLERR,"Insufficient memory on accelerator"); // communicate derivative of embedding function comm->forward_comm_pair(this); - + // compute forces on each atom on GPU - if (gpu_mode != GPU_FORCE) + if (gpu_mode != GPU_FORCE) eam_gpu_compute_force(NULL, eflag, vflag, eflag_atom, vflag_atom); else eam_gpu_compute_force(ilist, eflag, vflag, eflag_atom, vflag_atom); @@ -149,14 +149,14 @@ void PairEAMGPU::compute(int eflag, int vflag) void PairEAMGPU::init_style() { - if (force->newton_pair) + if (force->newton_pair) error->all(FLERR,"Cannot use newton pair with eam/gpu pair style"); - + // convert read-in file(s) to arrays and spline them file2array(); array2spline(); - + // Repeat cutsq calculation because done after call to init_style double maxcut = -1.0; double cut; @@ -173,18 +173,18 @@ void PairEAMGPU::init_style() } } double cell_size = sqrt(maxcut) + neighbor->skin; - + int maxspecial=0; if (atom->molecular) maxspecial=atom->maxspecial; int fp_size; int success = eam_gpu_init(atom->ntypes+1, cutforcesq, type2rhor, type2z2r, - type2frho, rhor_spline, z2r_spline, frho_spline, - rdr, rdrho, nrhor, nrho, nz2r, nfrho, nr, - atom->nlocal, atom->nlocal+atom->nghost, 300, - maxspecial, cell_size, gpu_mode, screen, fp_size); + type2frho, rhor_spline, z2r_spline, frho_spline, + rdr, rdrho, nrhor, nrho, nz2r, nfrho, nr, + atom->nlocal, atom->nlocal+atom->nghost, 300, + maxspecial, cell_size, gpu_mode, screen, fp_size); GPU_EXTRA::check_flag(success,error,world); - + if (gpu_mode == GPU_FORCE) { int irequest = neighbor->request(this); neighbor->requests[irequest]->half = 0; @@ -200,8 +200,8 @@ void PairEAMGPU::init_style() /* ---------------------------------------------------------------------- */ double PairEAMGPU::single(int i, int j, int itype, int jtype, - double rsq, double factor_coul, double factor_lj, - double &fforce) + double rsq, double factor_coul, double factor_lj, + double &fforce) { int m; double r,p,rhoip,rhojp,z2,z2p,recip,phi,phip,psip; @@ -213,7 +213,7 @@ double PairEAMGPU::single(int i, int j, int itype, int jtype, m = MIN(m,nr-1); p -= m; p = MIN(p,1.0); - + coeff = rhor_spline[type2rhor[itype][jtype]][m]; rhoip = (coeff[0]*p + coeff[1])*p + coeff[2]; coeff = rhor_spline[type2rhor[jtype][itype]][m]; @@ -221,7 +221,7 @@ double PairEAMGPU::single(int i, int j, int itype, int jtype, coeff = z2r_spline[type2z2r[itype][jtype]][m]; z2p = (coeff[0]*p + coeff[1])*p + coeff[2]; z2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - + double fp_i,fp_j; if (fp_single == false) { fp_i = ((double*)fp_pinned)[i]; @@ -230,20 +230,20 @@ double PairEAMGPU::single(int i, int j, int itype, int jtype, fp_i = ((float*)fp_pinned)[i]; fp_j = ((float*)fp_pinned)[j]; } - + recip = 1.0/r; phi = z2*recip; phip = z2p*recip - phi*recip; psip = fp_i*rhojp + fp_j*rhoip + phip; fforce = -psip*recip; - + return phi; } /* ---------------------------------------------------------------------- */ -int PairEAMGPU::pack_comm(int n, int *list, double *buf, int pbc_flag, - int *pbc) +int PairEAMGPU::pack_comm(int n, int *list, double *buf, int pbc_flag, + int *pbc) { int i,j,m; diff --git a/src/GPU/pair_eam_gpu.h b/src/GPU/pair_eam_gpu.h index 3b48c6b378..8766e5a5bb 100644 --- a/src/GPU/pair_eam_gpu.h +++ b/src/GPU/pair_eam_gpu.h @@ -5,7 +5,7 @@ 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 + 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. @@ -45,7 +45,7 @@ class PairEAMGPU : public PairEAM { double cpu_time; int *gpulist; void *fp_pinned; - bool fp_single; + bool fp_single; }; } diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index 9b8cca3628..9960df2c8f 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -378,7 +378,7 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype, if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]]; if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]]; } - + meff = mi*mj / (mi+mj); if (mask[i] & freeze_group_bit) meff = mj; if (mask[j] & freeze_group_bit) meff = mi; diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 6e3fddbe18..88075420b7 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -295,7 +295,7 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq, if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]]; if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]]; } - + meff = mi*mj / (mi+mj); if (mask[i] & freeze_group_bit) meff = mj; if (mask[j] & freeze_group_bit) meff = mi; diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index c7cd3ad3fc..1ac57ee2e7 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -587,7 +587,7 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype, double mi,mj,meff,damp,ccel,polyhertz; double vtr1,vtr2,vtr3,vrel,shrmag,rsht; double fs1,fs2,fs3,fs,fn; - + double *radius = atom->radius; radi = radius[i]; radj = radius[j]; @@ -658,7 +658,7 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype, if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]]; if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]]; } - + meff = mi*mj / (mi+mj); if (mask[i] & freeze_group_bit) meff = mj; if (mask[j] & freeze_group_bit) meff = mi; @@ -734,7 +734,7 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype, /* ---------------------------------------------------------------------- */ -int PairGranHookeHistory::pack_comm(int n, int *list, +int PairGranHookeHistory::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; diff --git a/src/KSPACE/pair_lj_cut_coul_long.cpp b/src/KSPACE/pair_lj_cut_coul_long.cpp index c5cfe503af..e83e50e239 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long.cpp @@ -1,1137 +1,1137 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Paul Crozier (SNL) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "pair_lj_cut_coul_long.h" -#include "atom.h" -#include "comm.h" -#include "force.h" -#include "kspace.h" -#include "update.h" -#include "integrate.h" -#include "respa.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "math_const.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace MathConst; - -#define EWALD_F 1.12837917 -#define EWALD_P 0.3275911 -#define A1 0.254829592 -#define A2 -0.284496736 -#define A3 1.421413741 -#define A4 -1.453152027 -#define A5 1.061405429 - -/* ---------------------------------------------------------------------- */ - -PairLJCutCoulLong::PairLJCutCoulLong(LAMMPS *lmp) : Pair(lmp) -{ - respa_enable = 1; - ftable = NULL; - qdist = 0.0; -} - -/* ---------------------------------------------------------------------- */ - -PairLJCutCoulLong::~PairLJCutCoulLong() -{ - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - - memory->destroy(cut_lj); - memory->destroy(cut_ljsq); - memory->destroy(epsilon); - memory->destroy(sigma); - memory->destroy(lj1); - memory->destroy(lj2); - memory->destroy(lj3); - memory->destroy(lj4); - memory->destroy(offset); - } - if (ftable) free_tables(); -} - -/* ---------------------------------------------------------------------- */ - -void PairLJCutCoulLong::compute(int eflag, int vflag) -{ - int i,ii,j,jj,inum,jnum,itype,jtype,itable; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double fraction,table; - double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; - double grij,expm2,prefactor,t,erfc; - int *ilist,*jlist,*numneigh,**firstneigh; - double rsq; - - evdwl = ecoul = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - double **x = atom->x; - double **f = atom->f; - double *q = atom->q; - int *type = atom->type; - int nlocal = atom->nlocal; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - qtmp = q[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_lj = special_lj[sbmask(j)]; - factor_coul = special_coul[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; - - if (rsq < cut_coulsq) { - if (!ncoultablebits || rsq <= tabinnersq) { - r = sqrt(rsq); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; - } else { - union_int_float_t rsq_lookup; - rsq_lookup.f = rsq; - itable = rsq_lookup.i & ncoulmask; - itable >>= ncoulshiftbits; - fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = qtmp*q[j] * table; - if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; - } - } - } else forcecoul = 0.0; - - if (rsq < cut_ljsq[itype][jtype]) { - r6inv = r2inv*r2inv*r2inv; - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); - } else forcelj = 0.0; - - fpair = (forcecoul + factor_lj*forcelj) * r2inv; - - 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) { - if (rsq < cut_coulsq) { - if (!ncoultablebits || rsq <= tabinnersq) - ecoul = prefactor*erfc; - else { - table = etable[itable] + fraction*detable[itable]; - ecoul = qtmp*q[j] * table; - } - if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; - } else ecoul = 0.0; - - if (rsq < cut_ljsq[itype][jtype]) { - evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - - offset[itype][jtype]; - evdwl *= factor_lj; - } else evdwl = 0.0; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); - } - } - } - - if (vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- */ - -void PairLJCutCoulLong::compute_inner() -{ - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; - double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; - - double **x = atom->x; - double **f = atom->f; - double *q = atom->q; - int *type = atom->type; - int nlocal = atom->nlocal; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; - - double cut_out_on = cut_respa[0]; - double cut_out_off = cut_respa[1]; - - double cut_out_diff = cut_out_off - cut_out_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - qtmp = q[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_lj = special_lj[sbmask(j)]; - factor_coul = special_coul[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - if (rsq < cut_out_off_sq) { - r2inv = 1.0/rsq; - forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; - - jtype = type[j]; - if (rsq < cut_ljsq[itype][jtype]) { - r6inv = r2inv*r2inv*r2inv; - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); - } else forcelj = 0.0; - - fpair = (forcecoul + factor_lj*forcelj) * r2inv; - if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.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; - } - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void PairLJCutCoulLong::compute_middle() -{ - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; - double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; - - double **x = atom->x; - double **f = atom->f; - double *q = atom->q; - int *type = atom->type; - int nlocal = atom->nlocal; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; - - double cut_in_off = cut_respa[0]; - double cut_in_on = cut_respa[1]; - double cut_out_on = cut_respa[2]; - double cut_out_off = cut_respa[3]; - - double cut_in_diff = cut_in_on - cut_in_off; - double cut_out_diff = cut_out_off - cut_out_on; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - qtmp = q[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_lj = special_lj[sbmask(j)]; - factor_coul = special_coul[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { - r2inv = 1.0/rsq; - forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; - - jtype = type[j]; - if (rsq < cut_ljsq[itype][jtype]) { - r6inv = r2inv*r2inv*r2inv; - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); - } else forcelj = 0.0; - - fpair = (forcecoul + factor_lj*forcelj) * r2inv; - if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); - } - if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.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; - } - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void PairLJCutCoulLong::compute_outer(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum,itype,jtype,itable; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double fraction,table; - double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; - double grij,expm2,prefactor,t,erfc; - double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; - double rsq; - - evdwl = ecoul = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = 0; - - double **x = atom->x; - double **f = atom->f; - double *q = atom->q; - int *type = atom->type; - int nlocal = atom->nlocal; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; - - double cut_in_off = cut_respa[2]; - double cut_in_on = cut_respa[3]; - - double cut_in_diff = cut_in_on - cut_in_off; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - qtmp = q[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_lj = special_lj[sbmask(j)]; - factor_coul = special_coul[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; - - if (rsq < cut_coulsq) { - if (!ncoultablebits || rsq <= tabinnersq) { - r = sqrt(rsq); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0); - if (rsq > cut_in_off_sq) { - if (rsq < cut_in_on_sq) { - rsw = (r - cut_in_off)/cut_in_diff; - forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw); - if (factor_coul < 1.0) - forcecoul -= - (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw); - } else { - forcecoul += prefactor; - if (factor_coul < 1.0) - forcecoul -= (1.0-factor_coul)*prefactor; - } - } - } else { - union_int_float_t rsq_lookup; - rsq_lookup.f = rsq; - itable = rsq_lookup.i & ncoulmask; - itable >>= ncoulshiftbits; - fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = qtmp*q[j] * table; - if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; - } - } - } else forcecoul = 0.0; - - if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) { - r6inv = r2inv*r2inv*r2inv; - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); - if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - forcelj *= rsw*rsw*(3.0 - 2.0*rsw); - } - } else forcelj = 0.0; - - fpair = (forcecoul + forcelj) * r2inv; - - 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) { - if (rsq < cut_coulsq) { - if (!ncoultablebits || rsq <= tabinnersq) { - ecoul = prefactor*erfc; - if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; - } else { - table = etable[itable] + fraction*detable[itable]; - ecoul = qtmp*q[j] * table; - if (factor_coul < 1.0) { - table = ptable[itable] + fraction*dptable[itable]; - prefactor = qtmp*q[j] * table; - ecoul -= (1.0-factor_coul)*prefactor; - } - } - } else ecoul = 0.0; - - if (rsq < cut_ljsq[itype][jtype]) { - r6inv = r2inv*r2inv*r2inv; - evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - - offset[itype][jtype]; - evdwl *= factor_lj; - } else evdwl = 0.0; - } - - if (vflag) { - if (rsq < cut_coulsq) { - if (!ncoultablebits || rsq <= tabinnersq) { - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; - } else { - table = vtable[itable] + fraction*dvtable[itable]; - forcecoul = qtmp*q[j] * table; - if (factor_coul < 1.0) { - table = ptable[itable] + fraction*dptable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; - } - } - } else forcecoul = 0.0; - - if (rsq <= cut_in_off_sq) { - r6inv = r2inv*r2inv*r2inv; - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); - } else if (rsq <= cut_in_on_sq) - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); - - fpair = (forcecoul + factor_lj*forcelj) * r2inv; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); - } - } - } -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::allocate() -{ - allocated = 1; - int n = atom->ntypes; - - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; - - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - - memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); - memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::settings(int narg, char **arg) -{ - if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command"); - - cut_lj_global = force->numeric(arg[0]); - if (narg == 1) cut_coul = cut_lj_global; - else cut_coul = force->numeric(arg[1]); - - // reset cutoffs that have been explicitly set - - if (allocated) { - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i+1; j <= atom->ntypes; j++) - if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; - } -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::coeff(int narg, char **arg) -{ - if (narg < 4 || narg > 5) - error->all(FLERR,"Incorrect args for pair coefficients"); - if (!allocated) allocate(); - - int ilo,ihi,jlo,jhi; - force->bounds(arg[0],atom->ntypes,ilo,ihi); - force->bounds(arg[1],atom->ntypes,jlo,jhi); - - double epsilon_one = force->numeric(arg[2]); - double sigma_one = force->numeric(arg[3]); - - double cut_lj_one = cut_lj_global; - if (narg == 5) cut_lj_one = force->numeric(arg[4]); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - epsilon[i][j] = epsilon_one; - sigma[i][j] = sigma_one; - cut_lj[i][j] = cut_lj_one; - setflag[i][j] = 1; - count++; - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::init_style() -{ - if (!atom->q_flag) - error->all(FLERR,"Pair style lj/cut/coul/long requires atom attribute q"); - - // request regular or rRESPA neighbor lists - - int irequest; - - if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { - int respa = 0; - if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; - if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; - - if (respa == 0) irequest = neighbor->request(this); - else if (respa == 1) { - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respaouter = 1; - } else { - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 2; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respamiddle = 1; - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respaouter = 1; - } - - } else irequest = neighbor->request(this); - - cut_coulsq = cut_coul * cut_coul; - - // set rRESPA cutoffs - - if (strstr(update->integrate_style,"respa") && - ((Respa *) update->integrate)->level_inner >= 0) - cut_respa = ((Respa *) update->integrate)->cutoff; - else cut_respa = NULL; - - // insure use of KSpace long-range solver, set g_ewald - - if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); - g_ewald = force->kspace->g_ewald; - - // setup force tables - - if (ncoultablebits) init_tables(); -} - -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairLJCutCoulLong::init_one(int i, int j) -{ - if (setflag[i][j] == 0) { - epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], - sigma[i][i],sigma[j][j]); - sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); - cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); - } - - // include TIP4P qdist in full cutoff, qdist = 0.0 if not TIP4P - - double cut = MAX(cut_lj[i][j],cut_coul+2.0*qdist); - cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; - - lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); - lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); - lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - - if (offset_flag) { - double ratio = sigma[i][j] / cut_lj[i][j]; - offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); - } else offset[i][j] = 0.0; - - cut_ljsq[j][i] = cut_ljsq[i][j]; - lj1[j][i] = lj1[i][j]; - lj2[j][i] = lj2[i][j]; - lj3[j][i] = lj3[i][j]; - lj4[j][i] = lj4[i][j]; - offset[j][i] = offset[i][j]; - - // check interior rRESPA cutoff - - if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) - error->all(FLERR,"Pair cutoff < Respa interior cutoff"); - - // compute I,J contribution to long-range tail correction - // count total # of atoms of type I and J via Allreduce - - if (tail_flag) { - int *type = atom->type; - int nlocal = atom->nlocal; - - double count[2],all[2]; - count[0] = count[1] = 0.0; - for (int k = 0; k < nlocal; k++) { - if (type[k] == i) count[0] += 1.0; - if (type[k] == j) count[1] += 1.0; - } - MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); - - double sig2 = sigma[i][j]*sigma[i][j]; - double sig6 = sig2*sig2*sig2; - double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; - double rc6 = rc3*rc3; - double rc9 = rc3*rc6; - etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); - ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); - } - - return cut; -} - -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::init_tables() -{ - int masklo,maskhi; - double r,grij,expm2,derfc,rsw; - double qqrd2e = force->qqrd2e; - - tabinnersq = tabinner*tabinner; - init_bitmap(tabinner,cut_coul,ncoultablebits, - masklo,maskhi,ncoulmask,ncoulshiftbits); - - int ntable = 1; - for (int i = 0; i < ncoultablebits; i++) ntable *= 2; - - // linear lookup tables of length N = 2^ncoultablebits - // stored value = value at lower edge of bin - // d values = delta from lower edge to upper edge of bin - - if (ftable) free_tables(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable*sizeof(double),"pair:vtable"); - memory->create(ptable,ntable*sizeof(double),"pair:ptable"); - memory->create(dvtable,ntable*sizeof(double),"pair:dvtable"); - memory->create(dptable,ntable*sizeof(double),"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - rsq_lookup.i = i << ncoulshiftbits; - rsq_lookup.i |= masklo; - if (rsq_lookup.f < tabinnersq) { - rsq_lookup.i = i << ncoulshiftbits; - rsq_lookup.i |= maskhi; - } - r = sqrtf(rsq_lookup.f); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - derfc = erfc(grij); - if (cut_respa == NULL) { - rtable[i] = rsq_lookup.f; - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - ctable[i] = qqrd2e/r; - etable[i] = qqrd2e/r * derfc; - } else { - rtable[i] = rsq_lookup.f; - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); - ctable[i] = 0.0; - etable[i] = qqrd2e/r * derfc; - ptable[i] = qqrd2e/r; - vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { - if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { - rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); - ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); - ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); - } else { - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - ctable[i] = qqrd2e/r; - } - } - } - minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - int ntablem1 = ntable - 1; - - for (int i = 0; i < ntablem1; i++) { - drtable[i] = 1.0/(rtable[i+1] - rtable[i]); - dftable[i] = ftable[i+1] - ftable[i]; - dctable[i] = ctable[i+1] - ctable[i]; - detable[i] = etable[i+1] - etable[i]; - } - - if (cut_respa) { - for (int i = 0; i < ntablem1; i++) { - dvtable[i] = vtable[i+1] - vtable[i]; - dptable[i] = ptable[i+1] - ptable[i]; - } - } - - // get the delta values for the last table entries - // tables are connected periodically between 0 and ntablem1 - - drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); - dftable[ntablem1] = ftable[0] - ftable[ntablem1]; - dctable[ntablem1] = ctable[0] - ctable[ntablem1]; - detable[ntablem1] = etable[0] - etable[ntablem1]; - if (cut_respa) { - dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; - dptable[ntablem1] = ptable[0] - ptable[ntablem1]; - } - - // get the correct delta values at itablemax - // smallest r is in bin itablemin - // largest r is in bin itablemax, which is itablemin-1, - // or ntablem1 if itablemin=0 - // deltas at itablemax only needed if corresponding rsq < cut*cut - // if so, compute deltas between rsq and cut*cut - - double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp; - itablemin = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - if (rsq_lookup.f < cut_coulsq) { - rsq_lookup.f = cut_coulsq; - r = sqrtf(rsq_lookup.f); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - derfc = erfc(grij); - - if (cut_respa == NULL) { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - e_tmp = qqrd2e/r * derfc; - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); - c_tmp = 0.0; - e_tmp = qqrd2e/r * derfc; - p_tmp = qqrd2e/r; - v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { - if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { - rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); - f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); - c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - } - } - } - - drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]); - dftable[itablemax] = f_tmp - ftable[itablemax]; - dctable[itablemax] = c_tmp - ctable[itablemax]; - detable[itablemax] = e_tmp - etable[itablemax]; - if (cut_respa) { - dvtable[itablemax] = v_tmp - vtable[itablemax]; - dptable[itablemax] = p_tmp - ptable[itablemax]; - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::write_restart(FILE *fp) -{ - write_restart_settings(fp); - - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); - if (setflag[i][j]) { - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut_lj[i][j],sizeof(double),1,fp); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::read_restart(FILE *fp) -{ - read_restart_settings(fp); - - allocate(); - - int i,j; - int me = comm->me; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) { - if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); - if (setflag[i][j]) { - if (me == 0) { - fread(&epsilon[i][j],sizeof(double),1,fp); - fread(&sigma[i][j],sizeof(double),1,fp); - fread(&cut_lj[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::write_restart_settings(FILE *fp) -{ - fwrite(&cut_lj_global,sizeof(double),1,fp); - fwrite(&cut_coul,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::read_restart_settings(FILE *fp) -{ - if (comm->me == 0) { - fread(&cut_lj_global,sizeof(double),1,fp); - fread(&cut_coul,sizeof(double),1,fp); - fread(&offset_flag,sizeof(int),1,fp); - fread(&mix_flag,sizeof(int),1,fp); - } - MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); -} - -/* ---------------------------------------------------------------------- - free memory for tables used in pair computations -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::free_tables() -{ - memory->destroy(rtable); - memory->destroy(drtable); - memory->destroy(ftable); - memory->destroy(dftable); - memory->destroy(ctable); - memory->destroy(dctable); - memory->destroy(etable); - memory->destroy(detable); - memory->destroy(vtable); - memory->destroy(dvtable); - memory->destroy(ptable); - memory->destroy(dptable); -} - -/* ---------------------------------------------------------------------- */ - -double PairLJCutCoulLong::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) -{ - double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor; - double fraction,table,forcecoul,forcelj,phicoul,philj; - int itable; - - r2inv = 1.0/rsq; - if (rsq < cut_coulsq) { - if (!ncoultablebits || rsq <= tabinnersq) { - r = sqrt(rsq); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; - } else { - union_int_float_t rsq_lookup_single; - rsq_lookup_single.f = rsq; - itable = rsq_lookup_single.i & ncoulmask; - itable >>= ncoulshiftbits; - fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = atom->q[i]*atom->q[j] * table; - if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = atom->q[i]*atom->q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; - } - } - } else forcecoul = 0.0; - - if (rsq < cut_ljsq[itype][jtype]) { - r6inv = r2inv*r2inv*r2inv; - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); - } else forcelj = 0.0; - - fforce = (forcecoul + factor_lj*forcelj) * r2inv; - - double eng = 0.0; - if (rsq < cut_coulsq) { - if (!ncoultablebits || rsq <= tabinnersq) - phicoul = prefactor*erfc; - else { - table = etable[itable] + fraction*detable[itable]; - phicoul = atom->q[i]*atom->q[j] * table; - } - if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; - eng += phicoul; - } - - if (rsq < cut_ljsq[itype][jtype]) { - philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - - offset[itype][jtype]; - eng += factor_lj*philj; - } - - return eng; -} - -/* ---------------------------------------------------------------------- */ - -void *PairLJCutCoulLong::extract(const char *str, int &dim) -{ - dim = 0; - if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; - dim = 2; - if (strcmp(str,"epsilon") == 0) return (void *) epsilon; - return NULL; -} +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lj_cut_coul_long.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "update.h" +#include "integrate.h" +#include "respa.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairLJCutCoulLong::PairLJCutCoulLong(LAMMPS *lmp) : Pair(lmp) +{ + respa_enable = 1; + ftable = NULL; + qdist = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +PairLJCutCoulLong::~PairLJCutCoulLong() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut_lj); + memory->destroy(cut_ljsq); + memory->destroy(epsilon); + memory->destroy(sigma); + memory->destroy(lj1); + memory->destroy(lj2); + memory->destroy(lj3); + memory->destroy(lj4); + memory->destroy(offset); + } + if (ftable) free_tables(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLong::compute(int eflag, int vflag) +{ + int i,ii,j,jj,inum,jnum,itype,jtype,itable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double fraction,table; + double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double grij,expm2,prefactor,t,erfc; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r2inv = 1.0/rsq; + + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + } else forcelj = 0.0; + + fpair = (forcecoul + factor_lj*forcelj) * r2inv; + + 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) { + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + ecoul = prefactor*erfc; + else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qtmp*q[j] * table; + } + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLong::compute_inner() +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; + double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double rsw; + int *ilist,*jlist,*numneigh,**firstneigh; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = listinner->inum; + ilist = listinner->ilist; + numneigh = listinner->numneigh; + firstneigh = listinner->firstneigh; + + double cut_out_on = cut_respa[0]; + double cut_out_off = cut_respa[1]; + + double cut_out_diff = cut_out_off - cut_out_on; + double cut_out_on_sq = cut_out_on*cut_out_on; + double cut_out_off_sq = cut_out_off*cut_out_off; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cut_out_off_sq) { + r2inv = 1.0/rsq; + forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; + + jtype = type[j]; + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + } else forcelj = 0.0; + + fpair = (forcecoul + factor_lj*forcelj) * r2inv; + if (rsq > cut_out_on_sq) { + rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; + fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.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; + } + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLong::compute_middle() +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; + double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double rsw; + int *ilist,*jlist,*numneigh,**firstneigh; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = listmiddle->inum; + ilist = listmiddle->ilist; + numneigh = listmiddle->numneigh; + firstneigh = listmiddle->firstneigh; + + double cut_in_off = cut_respa[0]; + double cut_in_on = cut_respa[1]; + double cut_out_on = cut_respa[2]; + double cut_out_off = cut_respa[3]; + + double cut_in_diff = cut_in_on - cut_in_off; + double cut_out_diff = cut_out_off - cut_out_on; + double cut_in_off_sq = cut_in_off*cut_in_off; + double cut_in_on_sq = cut_in_on*cut_in_on; + double cut_out_on_sq = cut_out_on*cut_out_on; + double cut_out_off_sq = cut_out_off*cut_out_off; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { + r2inv = 1.0/rsq; + forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; + + jtype = type[j]; + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + } else forcelj = 0.0; + + fpair = (forcecoul + factor_lj*forcelj) * r2inv; + if (rsq < cut_in_on_sq) { + rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; + fpair *= rsw*rsw*(3.0 - 2.0*rsw); + } + if (rsq > cut_out_on_sq) { + rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; + fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.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; + } + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLong::compute_outer(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype,itable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double fraction,table; + double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double grij,expm2,prefactor,t,erfc; + double rsw; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = listouter->inum; + ilist = listouter->ilist; + numneigh = listouter->numneigh; + firstneigh = listouter->firstneigh; + + double cut_in_off = cut_respa[2]; + double cut_in_on = cut_respa[3]; + + double cut_in_diff = cut_in_on - cut_in_off; + double cut_in_off_sq = cut_in_off*cut_in_off; + double cut_in_on_sq = cut_in_on*cut_in_on; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r2inv = 1.0/rsq; + + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0); + if (rsq > cut_in_off_sq) { + if (rsq < cut_in_on_sq) { + rsw = (r - cut_in_off)/cut_in_diff; + forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw); + if (factor_coul < 1.0) + forcecoul -= + (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw); + } else { + forcecoul += prefactor; + if (factor_coul < 1.0) + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq < cut_in_on_sq) { + rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; + forcelj *= rsw*rsw*(3.0 - 2.0*rsw); + } + } else forcelj = 0.0; + + fpair = (forcecoul + forcelj) * r2inv; + + 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) { + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + ecoul = prefactor*erfc; + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ptable[itable] + fraction*dptable[itable]; + prefactor = qtmp*q[j] * table; + ecoul -= (1.0-factor_coul)*prefactor; + } + } + } else ecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (vflag) { + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + table = vtable[itable] + fraction*dvtable[itable]; + forcecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ptable[itable] + fraction*dptable[itable]; + prefactor = qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq <= cut_in_off_sq) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + } else if (rsq <= cut_in_on_sq) + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + + fpair = (forcecoul + factor_lj*forcelj) * r2inv; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); + memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); + memory->create(epsilon,n+1,n+1,"pair:epsilon"); + memory->create(sigma,n+1,n+1,"pair:sigma"); + memory->create(lj1,n+1,n+1,"pair:lj1"); + memory->create(lj2,n+1,n+1,"pair:lj2"); + memory->create(lj3,n+1,n+1,"pair:lj3"); + memory->create(lj4,n+1,n+1,"pair:lj4"); + memory->create(offset,n+1,n+1,"pair:offset"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::settings(int narg, char **arg) +{ + if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command"); + + cut_lj_global = force->numeric(arg[0]); + if (narg == 1) cut_coul = cut_lj_global; + else cut_coul = force->numeric(arg[1]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::coeff(int narg, char **arg) +{ + if (narg < 4 || narg > 5) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double epsilon_one = force->numeric(arg[2]); + double sigma_one = force->numeric(arg[3]); + + double cut_lj_one = cut_lj_global; + if (narg == 5) cut_lj_one = force->numeric(arg[4]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + epsilon[i][j] = epsilon_one; + sigma[i][j] = sigma_one; + cut_lj[i][j] = cut_lj_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style lj/cut/coul/long requires atom attribute q"); + + // request regular or rRESPA neighbor lists + + int irequest; + + if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { + int respa = 0; + if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; + if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + + if (respa == 0) irequest = neighbor->request(this); + else if (respa == 1) { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } else { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 2; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respamiddle = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } + + } else irequest = neighbor->request(this); + + cut_coulsq = cut_coul * cut_coul; + + // set rRESPA cutoffs + + if (strstr(update->integrate_style,"respa") && + ((Respa *) update->integrate)->level_inner >= 0) + cut_respa = ((Respa *) update->integrate)->cutoff; + else cut_respa = NULL; + + // insure use of KSpace long-range solver, set g_ewald + + if (force->kspace == NULL) + error->all(FLERR,"Pair style is incompatible with KSpace style"); + g_ewald = force->kspace->g_ewald; + + // setup force tables + + if (ncoultablebits) init_tables(); +} + +/* ---------------------------------------------------------------------- + neighbor callback to inform pair style of neighbor list to use + regular or rRESPA +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::init_list(int id, NeighList *ptr) +{ + if (id == 0) list = ptr; + else if (id == 1) listinner = ptr; + else if (id == 2) listmiddle = ptr; + else if (id == 3) listouter = ptr; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLJCutCoulLong::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], + sigma[i][i],sigma[j][j]); + sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); + cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); + } + + // include TIP4P qdist in full cutoff, qdist = 0.0 if not TIP4P + + double cut = MAX(cut_lj[i][j],cut_coul+2.0*qdist); + cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; + + lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + + if (offset_flag) { + double ratio = sigma[i][j] / cut_lj[i][j]; + offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); + } else offset[i][j] = 0.0; + + cut_ljsq[j][i] = cut_ljsq[i][j]; + lj1[j][i] = lj1[i][j]; + lj2[j][i] = lj2[i][j]; + lj3[j][i] = lj3[i][j]; + lj4[j][i] = lj4[i][j]; + offset[j][i] = offset[i][j]; + + // check interior rRESPA cutoff + + if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) + error->all(FLERR,"Pair cutoff < Respa interior cutoff"); + + // compute I,J contribution to long-range tail correction + // count total # of atoms of type I and J via Allreduce + + if (tail_flag) { + int *type = atom->type; + int nlocal = atom->nlocal; + + double count[2],all[2]; + count[0] = count[1] = 0.0; + for (int k = 0; k < nlocal; k++) { + if (type[k] == i) count[0] += 1.0; + if (type[k] == j) count[1] += 1.0; + } + MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + + double sig2 = sigma[i][j]*sigma[i][j]; + double sig6 = sig2*sig2*sig2; + double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; + double rc6 = rc3*rc3; + double rc9 = rc3*rc6; + etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * + sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); + ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * + sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); + } + + return cut; +} + +/* ---------------------------------------------------------------------- + setup force tables used in compute routines +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::init_tables() +{ + int masklo,maskhi; + double r,grij,expm2,derfc,rsw; + double qqrd2e = force->qqrd2e; + + tabinnersq = tabinner*tabinner; + init_bitmap(tabinner,cut_coul,ncoultablebits, + masklo,maskhi,ncoulmask,ncoulshiftbits); + + int ntable = 1; + for (int i = 0; i < ncoultablebits; i++) ntable *= 2; + + // linear lookup tables of length N = 2^ncoultablebits + // stored value = value at lower edge of bin + // d values = delta from lower edge to upper edge of bin + + if (ftable) free_tables(); + + memory->create(rtable,ntable,"pair:rtable"); + memory->create(ftable,ntable,"pair:ftable"); + memory->create(ctable,ntable,"pair:ctable"); + memory->create(etable,ntable,"pair:etable"); + memory->create(drtable,ntable,"pair:drtable"); + memory->create(dftable,ntable,"pair:dftable"); + memory->create(dctable,ntable,"pair:dctable"); + memory->create(detable,ntable,"pair:detable"); + + if (cut_respa == NULL) { + vtable = ptable = dvtable = dptable = NULL; + } else { + memory->create(vtable,ntable*sizeof(double),"pair:vtable"); + memory->create(ptable,ntable*sizeof(double),"pair:ptable"); + memory->create(dvtable,ntable*sizeof(double),"pair:dvtable"); + memory->create(dptable,ntable*sizeof(double),"pair:dptable"); + } + + union_int_float_t rsq_lookup; + union_int_float_t minrsq_lookup; + int itablemin; + minrsq_lookup.i = 0 << ncoulshiftbits; + minrsq_lookup.i |= maskhi; + + for (int i = 0; i < ntable; i++) { + rsq_lookup.i = i << ncoulshiftbits; + rsq_lookup.i |= masklo; + if (rsq_lookup.f < tabinnersq) { + rsq_lookup.i = i << ncoulshiftbits; + rsq_lookup.i |= maskhi; + } + r = sqrtf(rsq_lookup.f); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + derfc = erfc(grij); + if (cut_respa == NULL) { + rtable[i] = rsq_lookup.f; + ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + ctable[i] = qqrd2e/r; + etable[i] = qqrd2e/r * derfc; + } else { + rtable[i] = rsq_lookup.f; + ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); + ctable[i] = 0.0; + etable[i] = qqrd2e/r * derfc; + ptable[i] = qqrd2e/r; + vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { + if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { + rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); + ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + } else { + ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + ctable[i] = qqrd2e/r; + } + } + } + minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); + } + + tabinnersq = minrsq_lookup.f; + + int ntablem1 = ntable - 1; + + for (int i = 0; i < ntablem1; i++) { + drtable[i] = 1.0/(rtable[i+1] - rtable[i]); + dftable[i] = ftable[i+1] - ftable[i]; + dctable[i] = ctable[i+1] - ctable[i]; + detable[i] = etable[i+1] - etable[i]; + } + + if (cut_respa) { + for (int i = 0; i < ntablem1; i++) { + dvtable[i] = vtable[i+1] - vtable[i]; + dptable[i] = ptable[i+1] - ptable[i]; + } + } + + // get the delta values for the last table entries + // tables are connected periodically between 0 and ntablem1 + + drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); + dftable[ntablem1] = ftable[0] - ftable[ntablem1]; + dctable[ntablem1] = ctable[0] - ctable[ntablem1]; + detable[ntablem1] = etable[0] - etable[ntablem1]; + if (cut_respa) { + dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; + dptable[ntablem1] = ptable[0] - ptable[ntablem1]; + } + + // get the correct delta values at itablemax + // smallest r is in bin itablemin + // largest r is in bin itablemax, which is itablemin-1, + // or ntablem1 if itablemin=0 + // deltas at itablemax only needed if corresponding rsq < cut*cut + // if so, compute deltas between rsq and cut*cut + + double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp; + itablemin = minrsq_lookup.i & ncoulmask; + itablemin >>= ncoulshiftbits; + int itablemax = itablemin - 1; + if (itablemin == 0) itablemax = ntablem1; + rsq_lookup.i = itablemax << ncoulshiftbits; + rsq_lookup.i |= maskhi; + + if (rsq_lookup.f < cut_coulsq) { + rsq_lookup.f = cut_coulsq; + r = sqrtf(rsq_lookup.f); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + derfc = erfc(grij); + + if (cut_respa == NULL) { + f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + c_tmp = qqrd2e/r; + e_tmp = qqrd2e/r * derfc; + } else { + f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); + c_tmp = 0.0; + e_tmp = qqrd2e/r * derfc; + p_tmp = qqrd2e/r; + v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { + if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { + rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); + f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + } else { + f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + c_tmp = qqrd2e/r; + } + } + } + + drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]); + dftable[itablemax] = f_tmp - ftable[itablemax]; + dctable[itablemax] = c_tmp - ctable[itablemax]; + detable[itablemax] = e_tmp - etable[itablemax]; + if (cut_respa) { + dvtable[itablemax] = v_tmp - vtable[itablemax]; + dptable[itablemax] = p_tmp - ptable[itablemax]; + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&epsilon[i][j],sizeof(double),1,fp); + fwrite(&sigma[i][j],sizeof(double),1,fp); + fwrite(&cut_lj[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::read_restart(FILE *fp) +{ + read_restart_settings(fp); + + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&epsilon[i][j],sizeof(double),1,fp); + fread(&sigma[i][j],sizeof(double),1,fp); + fread(&cut_lj[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::write_restart_settings(FILE *fp) +{ + fwrite(&cut_lj_global,sizeof(double),1,fp); + fwrite(&cut_coul,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_lj_global,sizeof(double),1,fp); + fread(&cut_coul,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- + free memory for tables used in pair computations +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::free_tables() +{ + memory->destroy(rtable); + memory->destroy(drtable); + memory->destroy(ftable); + memory->destroy(dftable); + memory->destroy(ctable); + memory->destroy(dctable); + memory->destroy(etable); + memory->destroy(detable); + memory->destroy(vtable); + memory->destroy(dvtable); + memory->destroy(ptable); + memory->destroy(dptable); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJCutCoulLong::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor; + double fraction,table,forcecoul,forcelj,phicoul,philj; + int itable; + + r2inv = 1.0/rsq; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup_single; + rsq_lookup_single.f = rsq; + itable = rsq_lookup_single.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = atom->q[i]*atom->q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = atom->q[i]*atom->q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + } else forcelj = 0.0; + + fforce = (forcecoul + factor_lj*forcelj) * r2inv; + + double eng = 0.0; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*erfc; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = atom->q[i]*atom->q[j] * table; + } + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + eng += phicoul; + } + + if (rsq < cut_ljsq[itype][jtype]) { + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + eng += factor_lj*philj; + } + + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJCutCoulLong::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + return NULL; +} diff --git a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp index d66ca1f5fc..3a3c2b11e1 100644 --- a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp @@ -614,4 +614,3 @@ double PairLJCutCoulLongTIP4P::memory_usage() bytes += 2 * nmax * sizeof(double); return bytes; } - diff --git a/src/MANYBODY/fix_qeq_comb.cpp b/src/MANYBODY/fix_qeq_comb.cpp index 364ae03cf9..6dbc12c521 100644 --- a/src/MANYBODY/fix_qeq_comb.cpp +++ b/src/MANYBODY/fix_qeq_comb.cpp @@ -5,7 +5,7 @@ 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 + 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. @@ -64,17 +64,17 @@ FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (strcmp(arg[iarg],"file") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/comb command"); if (me == 0) { - fp = fopen(arg[iarg+1],"w"); - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open fix qeq/comb file %s",arg[iarg+1]); - error->one(FLERR,str); - } + fp = fopen(arg[iarg+1],"w"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open fix qeq/comb file %s",arg[iarg+1]); + error->one(FLERR,str); + } } iarg += 2; } else error->all(FLERR,"Illegal fix qeq/comb command"); } - + nmax = atom->nmax; memory->create(qf,nmax,"qeq:qf"); memory->create(q1,nmax,"qeq:q1"); @@ -120,7 +120,7 @@ void FixQEQComb::init() error->all(FLERR,"Fix qeq/comb requires atom attribute q"); comb = (PairComb *) force->pair_match("comb",1); - if (comb == NULL) + if (comb == NULL) error->all(FLERR,"Must use pair_style comb with fix qeq/comb"); if (strstr(update->integrate_style,"respa")) @@ -176,10 +176,10 @@ void FixQEQComb::post_force(int vflag) memory->create(q2,nmax,"qeq:q2"); vector_atom = qf; } - + // more loops for first-time charge equilibrium - iloop = 0; + iloop = 0; if (firstflag) loopmax = 500; else loopmax = 200; @@ -187,15 +187,15 @@ void FixQEQComb::post_force(int vflag) if (me == 0 && fp) fprintf(fp,"Charge equilibration on step " BIGINT_FORMAT "\n", - update->ntimestep); - + update->ntimestep); + heatpq = 0.05; qmass = 0.016; dtq = 0.01; dtq2 = 0.5*dtq*dtq/qmass; double enegchk = 0.0; - double enegtot = 0.0; + double enegtot = 0.0; double enegmax = 0.0; double *q = atom->q; @@ -215,8 +215,8 @@ void FixQEQComb::post_force(int vflag) for (ii = 0; ii < inum; ii++) { i = ilist[ii]; if (mask[i] & groupbit) { - q1[i] += qf[i]*dtq2 - heatpq*q1[i]; - q[i] += q1[i]; + q1[i] += qf[i]*dtq2 - heatpq*q1[i]; + q[i] += q1[i]; } } comm->forward_comm_fix(this); @@ -224,14 +224,14 @@ void FixQEQComb::post_force(int vflag) if(comb) enegtot = comb->yasu_char(qf,igroup); enegtot /= ngroup; enegchk = enegmax = 0.0; - + for (ii = 0; ii < inum ; ii++) { i = ilist[ii]; if (mask[i] & groupbit) { - q2[i] = enegtot-qf[i]; - enegmax = MAX(enegmax,fabs(q2[i])); - enegchk += fabs(q2[i]); - qf[i] = q2[i]; + q2[i] = enegtot-qf[i]; + enegmax = MAX(enegmax,fabs(q2[i])); + enegchk += fabs(q2[i]); + qf[i] = q2[i]; } } @@ -239,27 +239,27 @@ void FixQEQComb::post_force(int vflag) enegchk = enegchkall/ngroup; MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world); enegmax = enegmaxall; - + if (enegchk <= precision && enegmax <= 100.0*precision) break; if (me == 0 && fp) fprintf(fp," iteration: %d, enegtot %.6g, " - "enegmax %.6g, fq deviation: %.6g\n", - iloop,enegtot,enegmax,enegchk); - + "enegmax %.6g, fq deviation: %.6g\n", + iloop,enegtot,enegmax,enegchk); + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; if (mask[i] & groupbit) - q1[i] += qf[i]*dtq2 - heatpq*q1[i]; + q1[i] += qf[i]*dtq2 - heatpq*q1[i]; } - } - + } + if (me == 0 && fp) { if (iloop == loopmax) fprintf(fp,"Charges did not converge in %d iterations\n",iloop); else fprintf(fp,"Charges converged in %d iterations to %.10f tolerance\n", - iloop,enegchk); + iloop,enegchk); } } @@ -291,18 +291,17 @@ int FixQEQComb::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) buf[m++] = atom->q[j]; } return 1; -} - -/* ---------------------------------------------------------------------- */ - -void FixQEQComb::unpack_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n ; - for (i = first; i < last; i++) atom->q[i] = buf[m++]; -} +} + +/* ---------------------------------------------------------------------- */ + +void FixQEQComb::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n ; + for (i = first; i < last; i++) atom->q[i] = buf[m++]; +} /* ---------------------------------------------------------------------- */ - diff --git a/src/MANYBODY/pair_bop.cpp b/src/MANYBODY/pair_bop.cpp index 2d027bc83c..a66fa34d41 100644 --- a/src/MANYBODY/pair_bop.cpp +++ b/src/MANYBODY/pair_bop.cpp @@ -19,7 +19,7 @@ The formulation for this work follows (a) D.G. Pettifor, et al., Mat. Sci. and Eng. A365, 2-13, (2004);(b) D.A. Murdick, et al., Phys. Rev. B 73, 045206 (2006);(c) D.G. Pettifor and I.I. Oleinik., Phys - Rev. Lett. 84, 4124 (2000); (d) D.K. Ward, et al., Phys. Rev. B 85, + Rev. Lett. 84, 4124 (2000); (d) D.K. Ward, et al., Phys. Rev. B 85, 115206 (2012). Copyright (2012) Sandia Corporation. Under the terms of Contract DE- @@ -289,7 +289,7 @@ void PairBOP::compute(int eflag, int vflag) int newton_pair = force->newton_pair; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - + inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -315,7 +315,7 @@ void PairBOP::compute(int eflag, int vflag) // Calculate Sigma Bond-Order - if(a_flag==1) { + if(a_flag==1) { if (otfly==0) sigmaBo_noa(); else sigmaBo_noa_otf(); } @@ -329,7 +329,7 @@ void PairBOP::compute(int eflag, int vflag) if (otfly==0) PiBo(); else PiBo_otf(); - n=0; + n=0; totE=0; for (ii = 0; ii < inum; ii++) { i=ilist[ii]; @@ -356,12 +356,12 @@ void PairBOP::compute(int eflag, int vflag) f[j][1]=f[j][1]-ftmp2; f[j][2]=f[j][2]-ftmp3; - // add repulsive and bond order components to total energy + // add repulsive and bond order components to total energy // (d) Eq.1 dE=-2.0*betaS[temp_ij]*sigB[n]-2.0*betaP[temp_ij]*piB[n]; totE+=dE+repul[temp_ij]; - if(evflag) { + if(evflag) { ev_tally_full(i,repul[temp_ij],dE,0.0,0.0,0.0,0.0); ev_tally_full(j,repul[temp_ij],dE,0.0,0.0,0.0,0.0); ev_tally_xyz(i,j,nlocal,newton_pair,0.0,0.0,-ftmp1,-ftmp2,-ftmp3, @@ -377,13 +377,13 @@ void PairBOP::compute(int eflag, int vflag) iij=itype*bop_types-itype*(itype+1)/2+jtype-1; else iij=jtype*bop_types-jtype*(jtype+1)/2+itype-1; - dis_ij[0]=x[j][0]-x[i][0]; - dis_ij[1]=x[j][1]-x[i][1]; - dis_ij[2]=x[j][2]-x[i][2]; + dis_ij[0]=x[j][0]-x[i][0]; + dis_ij[1]=x[j][1]-x[i][1]; + dis_ij[2]=x[j][2]-x[i][2]; rsq_ij=dis_ij[0]*dis_ij[0] +dis_ij[1]*dis_ij[1] +dis_ij[2]*dis_ij[2]; - r_ij=sqrt(rsq_ij); + r_ij=sqrt(rsq_ij); if(r_ij<=rcut[iij]) { ps=r_ij*rdr[iij]+1.0; ks=(int)ps; @@ -416,12 +416,12 @@ void PairBOP::compute(int eflag, int vflag) f[j][1]=f[j][1]-ftmp2; f[j][2]=f[j][2]-ftmp3; - // add repulsive and bond order components to total energy - // (d) Eq. 1 + // add repulsive and bond order components to total energy + // (d) Eq. 1 dE=-2.0*betaS_ij*sigB[n]-2.0*betaP_ij*piB[n]; totE+=dE+repul_ij; - if(evflag) { + if(evflag) { ev_tally_full(i,repul_ij,dE,0.0,0.0,0.0,0.0); ev_tally_full(j,repul_ij,dE,0.0,0.0,0.0,0.0); ev_tally_xyz(i,j,nlocal,newton_pair,0.0,0.0,-ftmp1,-ftmp2,-ftmp3, @@ -446,7 +446,7 @@ void PairBOP::allocate() { allocated = 1; int n = atom->ntypes; - + memory->create(rcut,npairs,"BOP:rcut"); memory->create(dr,npairs,"BOP:dr"); memory->create(rdr,npairs,"BOP:dr"); @@ -550,12 +550,12 @@ void PairBOP::coeff(int narg, char **arg) // match element names to BOP word types - if (me == 0) { + if (me == 0) { for (i = 3; i < narg; i++) { if (strcmp(arg[i],"NULL") == 0) { map[i-2] = -1; continue; - } + } for (j = 0; j < bop_types; j++) if (strcmp(arg[i],words[j]) == 0) break; map[i-2] = j; @@ -627,12 +627,12 @@ double PairBOP::init_one(int i, int j) int ii = map[i]+1; int jj = map[j]+1; - + int ij; if (ii==jj) ij=ii-1; else if (iiilist; firstneigh = list->firstneigh; numneigh = list->numneigh; - if(update_list!=0) + if(update_list!=0) memory_theta_grow(); else memory_theta_create(); @@ -724,7 +724,7 @@ void PairBOP::theta() else i=ii; itype = map[type[i]]+1; - + iilist=firstneigh[i]; for(jj=0;jj=npairs) { error->one(FLERR,"Too many atom pairs for pair bop"); } - disij[0][temp_ij]=x[j][0]-x[i][0]; - disij[1][temp_ij]=x[j][1]-x[i][1]; - disij[2][temp_ij]=x[j][2]-x[i][2]; + disij[0][temp_ij]=x[j][0]-x[i][0]; + disij[1][temp_ij]=x[j][1]-x[i][1]; + disij[2][temp_ij]=x[j][2]-x[i][2]; rsq=disij[0][temp_ij]*disij[0][temp_ij] +disij[1][temp_ij]*disij[1][temp_ij] +disij[2][temp_ij]*disij[2][temp_ij]; - rij[temp_ij]=sqrt(rsq); - if(rij[temp_ij]<=rcut[i12]) + rij[temp_ij]=sqrt(rsq); + if(rij[temp_ij]<=rcut[i12]) neigh_flag[temp_ij]=1; else neigh_flag[temp_ij]=0; @@ -789,7 +789,7 @@ void PairBOP::theta() error->one(FLERR,"Too many atom triplets for pair bop"); } temp_ik=BOP_index[i]+kk; - temp_ijk=cos_index[i]+n; + temp_ijk=cos_index[i]+n; if(temp_ijk>=cos_total) { error->one(FLERR,"Too many atom triplets for pair bop"); } @@ -817,7 +817,7 @@ void PairBOP::theta() dcAng[temp_ijk][2][1]=(disij[2][temp_ij]*rj1k1-cosAng[temp_ijk] *disij[2][temp_ik]*rj2)/(rj2k2); n++; - } + } } } } @@ -836,7 +836,7 @@ void PairBOP::theta_mod() /* The formulation differs slightly to avoid negative square roots in the calculation of Sigma^(1/2) of (a) Eq. 6 and (b) Eq. 11 */ - + void PairBOP::sigmaBo() { int nb_t,new_n_tot; @@ -895,7 +895,7 @@ void PairBOP::sigmaBo() int nall = nlocal+atom->nghost; firstneigh = list->firstneigh; numneigh = list->numneigh; - inum = list->inum; + inum = list->inum; ilist = list->ilist; n=0; @@ -913,7 +913,7 @@ void PairBOP::sigmaBo() create_sigma(nb_sg); for(itmp=0;itmptag; int newton_pair = force->newton_pair; int *type = atom->type; - + nlocal = atom->nlocal; int nall = nlocal+atom->nghost; firstneigh = list->firstneigh; numneigh = list->numneigh; - inum = list->inum; + inum = list->inum; ilist = list->ilist; n=0; - + //loop over all local atoms if(nb_sg>16) { @@ -2461,7 +2461,7 @@ void PairBOP::sigmaBo_noa() create_sigma(nb_sg); for(itmp=0;itmp=cos_total) { error->one(FLERR,"Too many atom triplets for pair bop"); } - gmean0=sigma_g0[itype-1][ktype-1][kptype-1]; - gmean1=sigma_g1[itype-1][ktype-1][kptype-1]; - gmean2=sigma_g2[itype-1][ktype-1][kptype-1]; + gmean0=sigma_g0[itype-1][ktype-1][kptype-1]; + gmean1=sigma_g1[itype-1][ktype-1][kptype-1]; + gmean2=sigma_g2[itype-1][ktype-1][kptype-1]; amean=cosAng[ang_ikkp]; gfactor2=gmean0+gmean1*amean +gmean2*amean*amean; @@ -2943,10 +2943,10 @@ void PairBOP::sigmaBo_noa() gfactor=gfactorsq*gfactorsq2; rfactorrt=betaS[temp_ik]*betaS[temp_kkp]; rfactor=rfactorrt*rfactorrt; - -//3rd CC is third term of Eq. 11 (c) for i atom + +//3rd CC is third term of Eq. 11 (c) for i atom //where j , k =neighbor of i & k' =neighbor of k - + CC=CC+gfactor*rfactor; } } @@ -2954,9 +2954,9 @@ void PairBOP::sigmaBo_noa() } } } - + //j is a neighbor of i and k is a neighbor of j not equal to i - + for(ktmp=0;ktmp-1)&&(bt_sg[m].j>-1)) { @@ -3235,7 +3235,7 @@ void PairBOP::sigmaBo_noa() +bt_sg[m].dBB[1]; bt_sg[m].dAAC[2]=bt_sg[m].dAA[2] +bt_sg[m].dBB[2]; - } + } } UT=EEC*FF+BBC+small3[iij]; UT=1.0/sqrt(UT); @@ -3243,7 +3243,7 @@ void PairBOP::sigmaBo_noa() // FFC is slightly modified form of (a) Eq. 31 // GGC is slightly modified form of (a) Eq. 32 // bndtmp is a slightly modified form of (a) Eq. 30 and (b) Eq. 8 - + bndtmp=(FF+sigma_delta[iij]*sigma_delta[iij]) +sigma_c[iij]*AAC+small4; UTcom=-0.5*UT*UT*UT; @@ -3284,7 +3284,7 @@ void PairBOP::sigmaBo_noa() } } } - + //This loop is to ensure there is not an error for atoms with no neighbors (deposition) if(nb_t==0) { @@ -3331,9 +3331,9 @@ void PairBOP::sigmaBo_noa() part2=1.0-part1*EE1/part0; part3=dsigB1*part1/part0; part4=part3/part0*EE1; - + // sigB is the final expression for (a) Eq. 6 and (b) Eq. 11 - + sigB[n]=dsigB1*part2; pp1=2.0*betaS[temp_ij]; for(m=0;mnlocal; int nall = nlocal + atom->nghost; - inum = list->inum; + inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; @@ -3465,12 +3465,12 @@ void PairBOP::sigmaBo_otf() if(allocate_sigma) { destroy_sigma(); } - + create_sigma(nb_sg); for(itmp=0;itmpx; double **f = atom->f; @@ -5377,7 +5377,7 @@ void PairBOP::sigmaBo_noa_otf() nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - inum = list->inum; + inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; @@ -5391,13 +5391,13 @@ void PairBOP::sigmaBo_noa_otf() } create_sigma(nb_sg); for(itmp=0;itmp-1)&&(bt_sg[m].j>-1)) { @@ -6399,15 +6399,15 @@ void PairBOP::sigmaBo_noa_otf() +bt_sg[m].dBB[1]; bt_sg[m].dAAC[2]=bt_sg[m].dAA[2] +bt_sg[m].dBB[2]; - } + } } UT=EEC*FF+BBC+small3[iij]; UT=1.0/sqrt(UT); - + // FFC is slightly modified form of (a) Eq. 31 // GGC is slightly modified form of (a) Eq. 32 // bndtmp is a slightly modified form of (a) Eq. 30 and (b) Eq. 8 - + FFC=BBC*UT; GGC=EEC*UT; bndtmp=(FF+sigma_delta[iij]*sigma_delta[iij]) @@ -6450,9 +6450,9 @@ void PairBOP::sigmaBo_noa_otf() } } } - + //This loop is to ensure there is not an error for atoms with no neighbors (deposition) - + if(nb_t==0) { if(j>i) { bt_sg[0].dSigB1[0]=bndtmp1*dis_ij[0]; @@ -6497,9 +6497,9 @@ void PairBOP::sigmaBo_noa_otf() part2=1.0-part1*EE1/part0; part3=dsigB1*part1/part0; part4=part3/part0*EE1; - + // sigB is the final expression for (a) Eq. 6 and (b) Eq. 11 - + sigB[n]=dsigB1*part2; pp1=2.0*betaS_ij; for(m=0;milist; numneigh = list->numneigh; firstneigh = list->firstneigh; - + n=0; // Loop over all local atoms for i @@ -6596,7 +6596,7 @@ void PairBOP::PiBo() if(allocate_pi) { destroy_pi(); } - create_pi(nb_pi); + create_pi(nb_pi); for(itmp=0;itmpnb_pi) { + if(nb_t>nb_pi) { new_n_tot=nb_pi+maxneigh; grow_pi(nb_pi,new_n_tot); nb_pi=new_n_tot; @@ -6696,7 +6696,7 @@ void PairBOP::PiBo() dbetaCapSq1=2.0*pi_p[itype-1]*betaS[temp_ik]*dBetaS[temp_ik] -2.0*betaP[temp_ik]*dBetaP[temp_ik]; -//AA is Eq. 37 (a) and Eq. 19 (b) or i atoms +//AA is Eq. 37 (a) and Eq. 19 (b) or i atoms //1st BB is first term of Eq. 38 (a) where j and k =neighbors i AA=AA+sinFactor*betaS[temp_ik]+cosFactor*betaP[temp_ik]; @@ -6745,7 +6745,7 @@ void PairBOP::PiBo() app2*dcAng[ang_jik][2][ngk] +agpdpr2*disij[2][temp_ik]; -// j and k and k' are different neighbors of i +// j and k and k' are different neighbors of i for(ltmp=0;ltmpntypes; + int n = atom->ntypes; nlocal = atom->nlocal; nghost = atom->nghost; nall = nlocal + nghost; @@ -9195,7 +9195,7 @@ double PairBOP::memory_usage() bytes += neigh_total * sizeof(double); // cos_index bytes += nall * sizeof(double); - } + } // pi_a bytes += npairs * sizeof(double); // pro_delta @@ -9240,9 +9240,9 @@ double PairBOP::memory_usage() bytes += npairs * sizeof(double); // r1 bytes += npairs * sizeof(double); -// sigma_beta0 +// sigma_beta0 bytes += npairs * sizeof(double); -// pi_beta0 +// pi_beta0 bytes += npairs * sizeof(double); // phi0 bytes += npairs * sizeof(double); @@ -9295,11 +9295,11 @@ double PairBOP::memory_usage() void PairBOP::memory_theta_create() { int nlocal,nghost,nall; - + nlocal = atom->nlocal; nghost = atom->nghost; nall = nlocal + nghost; - if(maxneigh<8) + if(maxneigh<8) neigh_ct=(maxneigh-1)*(maxneigh-1)*(maxneigh-1); else neigh_ct=(maxneigh-1)*(maxneigh-1); @@ -9333,11 +9333,11 @@ void PairBOP::memory_theta_create() void PairBOP::memory_theta_grow() { int nlocal,nghost,nall; - + nlocal = atom->nlocal; nghost = atom->nghost; nall = nlocal + nghost; - if(maxneigh<8) + if(maxneigh<8) neigh_ct=(maxneigh-1)*(maxneigh-1)*(maxneigh-1); else neigh_ct=(maxneigh-1)*(maxneigh-1); diff --git a/src/MANYBODY/pair_bop.h b/src/MANYBODY/pair_bop.h index c6a4557fc0..6fee56c292 100644 --- a/src/MANYBODY/pair_bop.h +++ b/src/MANYBODY/pair_bop.h @@ -5,12 +5,12 @@ 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 + certain rights in this software. This software is distributed under the GNU General Public License. The this work follows the formulation from (a) D.G. Pettifor, et al., Mat. Sci. and Eng. A365, 2-13, (2004) and (b) D.A. Murdick, et al., Phys. - Rev. B 73, 045206 (2006). (c) D.K. Ward, et al., Phys. Rev. B 85, 115206 + Rev. B 73, 045206 (2006). (c) D.K. Ward, et al., Phys. Rev. B 85, 115206 (2012) See the README file in the top-level LAMMPS directory. @@ -48,13 +48,13 @@ class PairBOP : public Pair { int maxbopn; // maximum size of bop neighbor list for allocation int maxnall; // maximum size of bop neighbor list for allocation int *map; // mapping from atom types to elements - int nelements; // # of unique elments + int nelements; // # of unique elments int nr; // increments for the BOP potential int nBOt; // second BO increments int bop_types; // number of elments in potential int npairs; // number of element pairs - char **elements; // names of unique elements - int ***elem2param; + char **elements; // names of unique elements + int ***elem2param; int nparams; int bop_step; int allocate_pi; @@ -69,7 +69,7 @@ class PairBOP : public Pair { int cos_total; // number of cosines stored if not using on the fly int neigh_ct; // limit for large arrays -/*Parameters variables*/ +/*Parameters variables*/ int ncutoff,nfunc; int a_flag; @@ -85,64 +85,64 @@ class PairBOP : public Pair { double which,alpha,alpha1,beta1,gamma1,alpha2,beta2,alpha3; double beta3,rsmall,rbig,rcore; char **words; - - double cutmax; //max cutoff for all elements - int otfly; //Defines whether to do on the fly + + double cutmax; //max cutoff for all elements + int otfly; //Defines whether to do on the fly //calculations of angles and distances //on the fly will slow down calculations //but requires less memory on = 1, off=0 - int table; //determines the method for reading in + int table; //determines the method for reading in //potential parameters a preset table //or generate the tables using a spline /* Neigh variables */ - + double *rcut,*dr,*rdr; double **disij,*rij; /*Triple variables */ - + double *cosAng,***dcosAng,***dcAng; - + /*Double variables */ - + double *betaS,*dBetaS,*betaP; double *dBetaP,*repul,*dRepul; - + /*Sigma variables */ - + int **itypeSigBk,*nSigBk; double *sigB; double *sigB1; - + /*Pi variables */ - + int **itypePiBk,*nPiBk; double *piB; - + /*Grids1 variables */ - + double **pBetaS,**pBetaS1,**pBetaS2,**pBetaS3; double **pBetaS4,**pBetaS5,**pBetaS6; - + /*Grids2 variables */ - + double **pBetaP,**pBetaP1,**pBetaP2,**pBetaP3; double **pBetaP4,**pBetaP5,**pBetaP6; - + /*Grids3 variables */ - + double **pRepul,**pRepul1,**pRepul2,**pRepul3; double **pRepul4,**pRepul5,**pRepul6; - + /*Grids4 variables */ - + double **FsigBO,**FsigBO1,**FsigBO2,**FsigBO3; double **FsigBO4,**FsigBO5,**FsigBO6; - double dBO,rdBO; - + double dBO,rdBO; + /* End of BOP variables */ double **rcmin,**rcmax,**rcmaxp; diff --git a/src/MOLECULE/atom_vec_angle.cpp b/src/MOLECULE/atom_vec_angle.cpp index ceac39e34c..8b1d2cd2f6 100644 --- a/src/MOLECULE/atom_vec_angle.cpp +++ b/src/MOLECULE/atom_vec_angle.cpp @@ -748,7 +748,7 @@ void AtomVecAngle::create_atom(int itype, double *coord) x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; mask[nlocal] = 1; - image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | + image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMAX; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; diff --git a/src/MOLECULE/atom_vec_bond.cpp b/src/MOLECULE/atom_vec_bond.cpp index 16a2577bda..375c0dc770 100644 --- a/src/MOLECULE/atom_vec_bond.cpp +++ b/src/MOLECULE/atom_vec_bond.cpp @@ -695,7 +695,7 @@ void AtomVecBond::create_atom(int itype, double *coord) x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; mask[nlocal] = 1; - image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | + image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMAX; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; diff --git a/src/MOLECULE/atom_vec_full.cpp b/src/MOLECULE/atom_vec_full.cpp index 2e180714b0..e7f009c569 100644 --- a/src/MOLECULE/atom_vec_full.cpp +++ b/src/MOLECULE/atom_vec_full.cpp @@ -894,7 +894,7 @@ void AtomVecFull::create_atom(int itype, double *coord) x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; mask[nlocal] = 1; - image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | + image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMAX; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; diff --git a/src/MOLECULE/atom_vec_molecular.cpp b/src/MOLECULE/atom_vec_molecular.cpp index 5d5cb0df21..bb21e5b0f4 100644 --- a/src/MOLECULE/atom_vec_molecular.cpp +++ b/src/MOLECULE/atom_vec_molecular.cpp @@ -878,7 +878,7 @@ void AtomVecMolecular::create_atom(int itype, double *coord) x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; mask[nlocal] = 1; - image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | + image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMAX; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; diff --git a/src/PERI/atom_vec_peri.cpp b/src/PERI/atom_vec_peri.cpp index 6c3992d774..8bccaf5dee 100644 --- a/src/PERI/atom_vec_peri.cpp +++ b/src/PERI/atom_vec_peri.cpp @@ -733,7 +733,7 @@ void AtomVecPeri::create_atom(int itype, double *coord) x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; mask[nlocal] = 1; - image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | + image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMAX; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; diff --git a/src/REPLICA/fix_event.cpp b/src/REPLICA/fix_event.cpp index c8d0273f24..f62f2983d8 100644 --- a/src/REPLICA/fix_event.cpp +++ b/src/REPLICA/fix_event.cpp @@ -108,7 +108,7 @@ void FixEvent::restore_event() // Since xevent is unwrapped coordinate, need to // adjust image flags when remapping - image[i] = ((tagint) IMGMAX << IMG2BITS) | + image[i] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMASK; domain->remap(x[i],image[i]); } diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp index 05170719ac..bd736607bc 100644 --- a/src/REPLICA/tad.cpp +++ b/src/REPLICA/tad.cpp @@ -839,7 +839,7 @@ void TAD::revert() x[i][0] = array_atom[i][0]; x[i][1] = array_atom[i][1]; x[i][2] = array_atom[i][2]; - image[i] = + image[i] = ((int(array_atom[i][5]) + (tagint) IMGMAX & IMGMASK) << IMG2BITS) | ((int(array_atom[i][4]) + (tagint) IMGMAX & IMGMASK) << IMGBITS) | (int(array_atom[i][3]) + IMGMAX & IMGMASK); diff --git a/src/USER-AWPMD/atom_vec_wavepacket.cpp b/src/USER-AWPMD/atom_vec_wavepacket.cpp index 538245d07a..e36741e59d 100644 --- a/src/USER-AWPMD/atom_vec_wavepacket.cpp +++ b/src/USER-AWPMD/atom_vec_wavepacket.cpp @@ -867,7 +867,7 @@ void AtomVecWavepacket::create_atom(int itype, double *coord) x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; mask[nlocal] = 1; - image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | + image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMAX; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; @@ -1001,7 +1001,7 @@ bigint AtomVecWavepacket::memory_usage() if (atom->memcheck("spin")) bytes += memory->usage(spin,nmax); if (atom->memcheck("eradius")) bytes += memory->usage(eradius,nmax); if (atom->memcheck("ervel")) bytes += memory->usage(ervel,nmax); - if (atom->memcheck("erforce")) + if (atom->memcheck("erforce")) bytes += memory->usage(erforce,nmax*comm->nthreads); if (atom->memcheck("ervelforce")) bytes += memory->usage(ervelforce,nmax); diff --git a/src/USER-EFF/atom_vec_electron.cpp b/src/USER-EFF/atom_vec_electron.cpp index ff1c472741..e1e3bc4104 100644 --- a/src/USER-EFF/atom_vec_electron.cpp +++ b/src/USER-EFF/atom_vec_electron.cpp @@ -731,7 +731,7 @@ void AtomVecElectron::create_atom(int itype, double *coord) x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; mask[nlocal] = 1; - image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | + image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMAX; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; diff --git a/src/USER-MISC/bond_harmonic_shift.cpp b/src/USER-MISC/bond_harmonic_shift.cpp index 21b0ff40ab..3c2a060831 100644 --- a/src/USER-MISC/bond_harmonic_shift.cpp +++ b/src/USER-MISC/bond_harmonic_shift.cpp @@ -83,7 +83,7 @@ void BondHarmonicShift::compute(int eflag, int vflag) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - if (eflag) + if (eflag) ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) ); // apply force to each of 2 atoms diff --git a/src/USER-MISC/bond_harmonic_shift_cut.cpp b/src/USER-MISC/bond_harmonic_shift_cut.cpp index d867b5e0cf..bcaac0990f 100644 --- a/src/USER-MISC/bond_harmonic_shift_cut.cpp +++ b/src/USER-MISC/bond_harmonic_shift_cut.cpp @@ -85,7 +85,7 @@ void BondHarmonicShiftCut::compute(int eflag, int vflag) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - if (eflag) + if (eflag) ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type])); // apply force to each of 2 atoms diff --git a/src/USER-MOLFILE/dump_molfile.cpp b/src/USER-MOLFILE/dump_molfile.cpp index 1d01884302..90d15f6779 100644 --- a/src/USER-MOLFILE/dump_molfile.cpp +++ b/src/USER-MOLFILE/dump_molfile.cpp @@ -5,7 +5,7 @@ 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 + 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. @@ -41,11 +41,11 @@ typedef MolfileInterface MFI; // path // template (import name and topology information from file) // bonds (write out bond information) -// topology (write out all topology information) +// topology (write out all topology information) /* ---------------------------------------------------------------------- */ -DumpMolfile::DumpMolfile(LAMMPS *lmp, int narg, char **arg) +DumpMolfile::DumpMolfile(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) { if (narg < 6) error->all(FLERR,"Illegal dump molfile command"); @@ -81,9 +81,9 @@ DumpMolfile::DumpMolfile(LAMMPS *lmp, int narg, char **arg) // allocate global array for atom coords bigint n = group->count(igroup); - if (n > MAXSMALLINT/sizeof(float)) + if (n > MAXSMALLINT/sizeof(float)) error->all(FLERR,"Too many atoms for dump molfile"); - if (n < 1) + if (n < 1) error->all(FLERR,"Not enough atoms for dump molfile"); natoms = static_cast(n); @@ -127,7 +127,7 @@ DumpMolfile::~DumpMolfile() memory->destroy(radiuses); delete mf; } - + if (typenames) { for (int i = 1; i <= ntypes; i++) delete [] typenames[i]; @@ -146,7 +146,7 @@ void DumpMolfile::init_style() if (me == 0) { - /* initialize typenames array to numeric types by default */ + /* initialize typenames array to numeric types by default */ if (typenames == NULL) { typenames = new char*[ntypes+1]; for (int itype = 1; itype <= ntypes; itype++) { @@ -249,7 +249,7 @@ void DumpMolfile::write() write_data(nlines,buf); } - + } else { MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status); MPI_Rsend(buf,nme*size_one,MPI_DOUBLE,0,0,world); @@ -282,7 +282,7 @@ void DumpMolfile::openfile() char *p2 = filecurrent; while (p1 != ptr) *p2++ = *p1++; - + if (padflag == 0) { sprintf(p2,BIGINT_FORMAT "%s",update->ntimestep,ptr+1); } else { @@ -393,7 +393,7 @@ void DumpMolfile::write_data(int n, double *mybuf) if (atom->rmass_flag) { mf->property(MFI::P_MASS,masses); - } else { + } else { mf->property(MFI::P_MASS,types,atom->mass); } diff --git a/src/USER-MOLFILE/dump_molfile.h b/src/USER-MOLFILE/dump_molfile.h index 964ff349bc..d14b4f58c4 100644 --- a/src/USER-MOLFILE/dump_molfile.h +++ b/src/USER-MOLFILE/dump_molfile.h @@ -5,7 +5,7 @@ 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 + 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. diff --git a/src/USER-MOLFILE/molfile_interface.cpp b/src/USER-MOLFILE/molfile_interface.cpp index 13bb49fd57..f3b63e7eda 100644 --- a/src/USER-MOLFILE/molfile_interface.cpp +++ b/src/USER-MOLFILE/molfile_interface.cpp @@ -5,7 +5,7 @@ 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 + 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. @@ -50,7 +50,7 @@ extern "C" { } plugin_reginfo_t; // callback function for plugin registration. - static int plugin_register_cb(void *v, vmdplugin_t *p) + static int plugin_register_cb(void *v, vmdplugin_t *p) { plugin_reginfo_t *r = static_cast(v); // make sure we have the proper plugin type (native reader) @@ -63,10 +63,10 @@ extern "C" { } /* periodic table of elements for translation of ordinal to atom type */ - static const char *pte_label[] = { + static const char *pte_label[] = { "X", "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P" , "S", "Cl", "Ar", "K", "Ca", "Sc", - "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", + "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", @@ -79,22 +79,22 @@ extern "C" { static const int nr_pte_entries = sizeof(pte_label) / sizeof(char *); /* corresponding table of masses. */ - static const float pte_mass[] = { - /* X */ 0.00000, 1.00794, 4.00260, 6.941, 9.012182, 10.811, - /* C */ 12.0107, 14.0067, 15.9994, 18.9984032, 20.1797, + static const float pte_mass[] = { + /* X */ 0.00000, 1.00794, 4.00260, 6.941, 9.012182, 10.811, + /* C */ 12.0107, 14.0067, 15.9994, 18.9984032, 20.1797, /* Na */ 22.989770, 24.3050, 26.981538, 28.0855, 30.973761, /* S */ 32.065, 35.453, 39.948, 39.0983, 40.078, 44.955910, /* Ti */ 47.867, 50.9415, 51.9961, 54.938049, 55.845, 58.9332, - /* Ni */ 58.6934, 63.546, 65.409, 69.723, 72.64, 74.92160, - /* Se */ 78.96, 79.904, 83.798, 85.4678, 87.62, 88.90585, + /* Ni */ 58.6934, 63.546, 65.409, 69.723, 72.64, 74.92160, + /* Se */ 78.96, 79.904, 83.798, 85.4678, 87.62, 88.90585, /* Zr */ 91.224, 92.90638, 95.94, 98.0, 101.07, 102.90550, - /* Pd */ 106.42, 107.8682, 112.411, 114.818, 118.710, 121.760, - /* Te */ 127.60, 126.90447, 131.293, 132.90545, 137.327, + /* Pd */ 106.42, 107.8682, 112.411, 114.818, 118.710, 121.760, + /* Te */ 127.60, 126.90447, 131.293, 132.90545, 137.327, /* La */ 138.9055, 140.116, 140.90765, 144.24, 145.0, 150.36, - /* Eu */ 151.964, 157.25, 158.92534, 162.500, 164.93032, + /* Eu */ 151.964, 157.25, 158.92534, 162.500, 164.93032, /* Er */ 167.259, 168.93421, 173.04, 174.967, 178.49, 180.9479, - /* W */ 183.84, 186.207, 190.23, 192.217, 195.078, 196.96655, - /* Hg */ 200.59, 204.3833, 207.2, 208.98038, 209.0, 210.0, 222.0, + /* W */ 183.84, 186.207, 190.23, 192.217, 195.078, 196.96655, + /* Hg */ 200.59, 204.3833, 207.2, 208.98038, 209.0, 210.0, 222.0, /* Fr */ 223.0, 226.0, 227.0, 232.0381, 231.03588, 238.02891, /* Np */ 237.0, 244.0, 243.0, 247.0, 247.0, 251.0, 252.0, 257.0, /* Md */ 258.0, 259.0, 262.0, 261.0, 262.0, 266.0, 264.0, 269.0, @@ -103,24 +103,24 @@ extern "C" { /* * corresponding table of VDW radii. - * van der Waals radii are taken from A. Bondi, - * J. Phys. Chem., 68, 441 - 452, 1964, - * except the value for H, which is taken from R.S. Rowland & R. Taylor, - * J.Phys.Chem., 100, 7384 - 7391, 1996. Radii that are not available in + * van der Waals radii are taken from A. Bondi, + * J. Phys. Chem., 68, 441 - 452, 1964, + * except the value for H, which is taken from R.S. Rowland & R. Taylor, + * J.Phys.Chem., 100, 7384 - 7391, 1996. Radii that are not available in * either of these publications have RvdW = 2.00 Å. - * The radii for Ions (Na, K, Cl, Ca, Mg, and Cs are based on the CHARMM27 + * The radii for Ions (Na, K, Cl, Ca, Mg, and Cs are based on the CHARMM27 * Rmin/2 parameters for (SOD, POT, CLA, CAL, MG, CES) by default. */ - static const float pte_vdw_radius[] = { - /* X */ 1.5, 1.2, 1.4, 1.82, 2.0, 2.0, - /* C */ 1.7, 1.55, 1.52, 1.47, 1.54, + static const float pte_vdw_radius[] = { + /* X */ 1.5, 1.2, 1.4, 1.82, 2.0, 2.0, + /* C */ 1.7, 1.55, 1.52, 1.47, 1.54, /* Na */ 1.36, 1.18, 2.0, 2.1, 1.8, /* S */ 1.8, 2.27, 1.88, 1.76, 1.37, 2.0, /* Ti */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, /* Ni */ 1.63, 1.4, 1.39, 1.07, 2.0, 1.85, - /* Se */ 1.9, 1.85, 2.02, 2.0, 2.0, 2.0, + /* Se */ 1.9, 1.85, 2.02, 2.0, 2.0, 2.0, /* Zr */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, - /* Pd */ 1.63, 1.72, 1.58, 1.93, 2.17, 2.0, + /* Pd */ 1.63, 1.72, 1.58, 1.93, 2.17, 2.0, /* Te */ 2.06, 1.98, 2.16, 2.1, 2.0, /* La */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, /* Eu */ 2.0, 2.0, 2.0, 2.0, 2.0, @@ -166,7 +166,7 @@ extern "C" { case 8: return 1.3; /* O */ case 9: return 1.2; /* F */ case 15: return 1.5; /* P */ - case 16: return 1.9; /* S */ + case 16: return 1.9; /* S */ } #endif @@ -177,13 +177,13 @@ extern "C" { { int i; char atom[3]; - + /* zap string */ atom[0] = (char) 0; atom[1] = (char) 0; atom[2] = (char) 0; - /* if we don't have a null-pointer, there must be at least two - * chars, which is all we need. we convert to the capitalization + /* if we don't have a null-pointer, there must be at least two + * chars, which is all we need. we convert to the capitalization * convention of the table above during assignment. */ if (label != NULL) { atom[0] = (char) toupper((int) label[0]); @@ -191,12 +191,12 @@ extern "C" { } /* discard numbers in atom label */ if (isdigit(atom[1])) atom[1] = (char) 0; - + for (i=0; i < nr_pte_entries; ++i) { if ( (pte_label[i][0] == atom[0]) && (pte_label[i][1] == atom[1]) ) return i; } - + return 0; } @@ -217,12 +217,12 @@ extern "C" { if (ind < 1) return 0; /* no non-whitespace characters */ - + for (i=0; i < nr_pte_entries; ++i) { - if ((toupper(pte_label[i][0]) == atom[0]) && (toupper(pte_label[i][1]) == atom[1])) + if ((toupper(pte_label[i][0]) == atom[0]) && (toupper(pte_label[i][1]) == atom[1])) return i; } - } + } return 0; } @@ -246,7 +246,7 @@ extern "C" { dirhandle_t *d; int len; - if (dirname == NULL) + if (dirname == NULL) return NULL; d = new dirhandle_t; @@ -275,9 +275,9 @@ extern "C" { static char *my_readdir(dirhandle_t *d) { if (FindNextFile(d->h, &(d->fd))) { - return d->fd.cFileName; + return d->fd.cFileName; } - return NULL; + return NULL; } // close directory handle @@ -298,10 +298,10 @@ extern "C" { // report error message from dlopen static const char *my_dlerror(void) { - static CHAR szBuf[80]; + static CHAR szBuf[80]; DWORD dw = GetLastError(); - sprintf(szBuf, "my_dlopen failed: GetLastError returned %u\n", dw); + sprintf(szBuf, "my_dlopen failed: GetLastError returned %u\n", dw); return szBuf; } @@ -645,7 +645,7 @@ int MolfileInterface::open(const char *name, int *natoms) if (!_plugin || !_dso || !natoms) return E_FILE; molfile_plugin_t *p = static_cast(_plugin); - + if (_mode & M_WRITE) _ptr = p->open_file_write(name,_type,*natoms); else if (_mode & M_READ) @@ -688,7 +688,7 @@ int MolfileInterface::structure() optflags |= (_props & P_CHRG) ? MOLFILE_CHARGE : 0; optflags |= (_props & P_RADS) ? MOLFILE_RADIUS : 0; optflags |= (_props & P_ATMN) ? MOLFILE_ATOMICNUMBER : 0; - + molfile_atom_t *a = static_cast(_info); p->write_structure(_ptr,optflags,a); } else if (_mode & M_RSTRUCT) { @@ -743,7 +743,7 @@ int MolfileInterface::timestep(float *coords, float *vels, molfile_plugin_t *p = static_cast(_plugin); molfile_timestep_t *t = new molfile_timestep_t; int rv; - + if (_mode & M_WRITE) { t->coords = coords; t->velocities = vels; @@ -845,7 +845,7 @@ static int read_int_property(molfile_atom_t &a, const int propid) PROPUPDATE(MolfileInterface::P_RESI,resid,prop); PROPUPDATE(MolfileInterface::P_ATMN,atomicnumber,prop); - + PROPUPDATE((MolfileInterface::P_ATMN|MolfileInterface::P_NAME), name,sprop); PROPUPDATE((MolfileInterface::P_ATMN|MolfileInterface::P_TYPE), @@ -896,7 +896,7 @@ static const char *read_string_property(molfile_atom_t &a, // floating point props static int write_atom_property(molfile_atom_t &a, const int propid, - const float prop) + const float prop) { int plist = MolfileInterface::P_NONE; PROPUPDATE(MolfileInterface::P_OCCP,occupancy,prop); @@ -910,7 +910,7 @@ static int write_atom_property(molfile_atom_t &a, // double precision floating point props static int write_atom_property(molfile_atom_t &a, const int propid, - const double prop) + const double prop) { return write_atom_property(a,propid,static_cast(prop)); } @@ -918,7 +918,7 @@ static int write_atom_property(molfile_atom_t &a, // integer and derived props static int write_atom_property(molfile_atom_t &a, const int propid, - const int prop) + const int prop) { int plist = MolfileInterface::P_NONE; PROPUPDATE(MolfileInterface::P_RESI,resid,prop); @@ -935,7 +935,7 @@ static int write_atom_property(molfile_atom_t &a, // integer and derived props static int write_atom_property(molfile_atom_t &a, const int propid, - const char *prop) + const char *prop) { int plist = MolfileInterface::P_NONE; PROPSTRCPY(MolfileInterface::P_NAME,name,prop); @@ -1078,13 +1078,13 @@ int MolfileInterface::property(int propid, double *prop) _props |= write_atom_property(a[IDX],P_CHAI,buf) // set/get atom integer property -int MolfileInterface::property(int propid, int idx, int *prop) +int MolfileInterface::property(int propid, int idx, int *prop) { if ((_info == NULL) || (prop == NULL) || (idx < 0) || (idx >= _natoms)) return P_NONE; molfile_atom_t *a = static_cast(_info); - + if (_mode & M_WSTRUCT) { char buf[64]; @@ -1103,13 +1103,13 @@ int MolfileInterface::property(int propid, int idx, int *prop) } // set/get per type integer property -int MolfileInterface::property(int propid, int *types, int *prop) +int MolfileInterface::property(int propid, int *types, int *prop) { if ((_info == NULL) || (types == NULL) || (prop == NULL)) return P_NONE; molfile_atom_t *a = static_cast(_info); - + if (_mode & M_WSTRUCT) { char buf[64]; @@ -1132,13 +1132,13 @@ int MolfileInterface::property(int propid, int *types, int *prop) } // set/get per atom integer property -int MolfileInterface::property(int propid, int *prop) +int MolfileInterface::property(int propid, int *prop) { if ((_info == NULL) || (prop == NULL)) return P_NONE; molfile_atom_t *a = static_cast(_info); - + if (_mode & M_WSTRUCT) { char buf[64]; @@ -1150,7 +1150,7 @@ int MolfileInterface::property(int propid, int *prop) sprintf(buf,"%d",prop[i]); INT_TO_STRING_BODY(i); } - } + } } if (_mode & M_RSTRUCT) { @@ -1169,7 +1169,7 @@ int MolfileInterface::property(int propid, int idx, char *prop) return P_NONE; molfile_atom_t *a = static_cast(_info); - + if (_mode & M_WSTRUCT) { _props |= write_atom_property(a[idx], propid, prop); } @@ -1187,7 +1187,7 @@ int MolfileInterface::property(int propid, int *types, char **prop) return P_NONE; molfile_atom_t *a = static_cast(_info); - + if (_mode & M_WSTRUCT) { for (int i=0; i < _natoms; ++i) { _props |= write_atom_property(a[i], propid, prop[types[i]]); @@ -1208,7 +1208,7 @@ int MolfileInterface::property(int propid, char **prop) return P_NONE; molfile_atom_t *a = static_cast(_info); - + if (_mode & M_WSTRUCT) { for (int i=0; i < _natoms; ++i) { _props |= write_atom_property(a[i], propid, prop[i]); @@ -1221,5 +1221,3 @@ int MolfileInterface::property(int propid, char **prop) return _props; } - - diff --git a/src/USER-MOLFILE/molfile_interface.h b/src/USER-MOLFILE/molfile_interface.h index f7dc4f5221..1146fb8251 100644 --- a/src/USER-MOLFILE/molfile_interface.h +++ b/src/USER-MOLFILE/molfile_interface.h @@ -5,7 +5,7 @@ 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 + 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. @@ -23,7 +23,7 @@ namespace LAMMPS_NS { // This class provides an abstract interface // to the VMD plugin library. -class MolfileInterface +class MolfileInterface { public: // plugin modes. diff --git a/src/USER-MOLFILE/molfile_plugin.h b/src/USER-MOLFILE/molfile_plugin.h index 94f6dc36ac..7a2d7ca42e 100644 --- a/src/USER-MOLFILE/molfile_plugin.h +++ b/src/USER-MOLFILE/molfile_plugin.h @@ -15,10 +15,10 @@ * ***************************************************************************/ -/** @file +/** @file * API for C extensions to define a way to load structure, coordinate, - * trajectory, and volumetric data files - */ + * trajectory, and volumetric data files + */ #ifndef MOL_FILE_PLUGIN_H #define MOL_FILE_PLUGIN_H @@ -74,16 +74,16 @@ typedef struct { } molfile_metadata_t; -/* - * Struct for specifying atoms in a molecular structure. The first - * six components are required, the rest are optional and their presence is +/* + * Struct for specifying atoms in a molecular structure. The first + * six components are required, the rest are optional and their presence is * indicating by setting the corresponding bit in optsflag. When omitted, - * the application (for read_structure) or plugin (for write_structure) - * must be able to supply default values if the missing parameters are + * the application (for read_structure) or plugin (for write_structure) + * must be able to supply default values if the missing parameters are * part of its internal data structure. * Note that it is not possible to specify coordinates with this structure. - * This is intentional; all coordinate I/O is done with the read_timestep and - * write_timestep functions. + * This is intentional; all coordinate I/O is done with the read_timestep and + * write_timestep functions. */ /** @@ -128,7 +128,7 @@ typedef struct { #define MOLFILE_CTNUMBER 0x0200 /**< ctnumber provided */ #endif #define MOLFILE_BADOPTIONS 0xFFFFFFFF /**< Detect badly behaved plugins */ - + /*@}*/ /*@{*/ @@ -150,7 +150,7 @@ typedef struct molfile_timestep_metadata { /* * Per-timestep atom coordinates and periodic cell information - */ + */ typedef struct { float *coords; /**< coordinates of all atoms, arranged xyzxyzxyz */ #if vmdplugin_ABIVERSION > 10 @@ -158,16 +158,16 @@ typedef struct { /**< NULL unless has_velocities is set */ #endif - /*@{*/ + /*@{*/ /** * Unit cell specification of the form A, B, C, alpha, beta, gamma. * notes: A, B, C are side lengths of the unit cell * alpha = angle between b and c * beta = angle between a and c * gamma = angle between a and b - */ - float A, B, C, alpha, beta, gamma; - /*@}*/ + */ + float A, B, C, alpha, beta, gamma; + /*@}*/ #if vmdplugin_ABIVERSION > 10 double physical_time; /**< physical time point associated with this frame */ @@ -188,7 +188,7 @@ typedef struct { /** * Metadata for volumetric datasets, read initially and used for subsequent - * memory allocations and file loading. + * memory allocations and file loading. */ typedef struct { char dataname[256]; /**< name of volumetric data set */ @@ -198,20 +198,20 @@ typedef struct { * x/y/z axis: * These the three cell sides, providing both direction and length * (not unit vectors) for the x, y, and z axes. In the simplest - * case, these would be <0,size,0> and <0,0,size) for + * case, these would be <0,size,0> and <0,0,size) for * an orthogonal cubic volume set. For other cell shapes these * axes can be oriented non-orthogonally, and the parallelpiped * may have different side lengths, not just a cube/rhombus. */ - float xaxis[3]; /**< direction (and length) for X axis */ + float xaxis[3]; /**< direction (and length) for X axis */ float yaxis[3]; /**< direction (and length) for Y axis */ float zaxis[3]; /**< direction (and length) for Z axis */ /* - * x/y/z size: + * x/y/z size: * Number of grid cells along each axis. This is _not_ the * physical size of the box, this is the number of voxels in each - * direction, independent of the shape of the volume set. + * direction, independent of the shape of the volume set. */ int xsize; /**< number of grid cells along the X axis */ int ysize; /**< number of grid cells along the Y axis */ @@ -297,10 +297,10 @@ typedef struct { /** * QM run info. Parameters that stay unchanged during a single file. - */ + */ typedef struct { int nproc; /**< number of processors used. */ - int memory; /**< amount of memory used in Mbyte. */ + int memory; /**< amount of memory used in Mbyte. */ int runtype; /**< flag indicating the calculation method. */ int scftype; /**< SCF type: RHF, UHF, ROHF, GVB or MCSCF wfn. */ int status; /**< indicates wether SCF and geometry optimization @@ -333,9 +333,9 @@ typedef struct { * array size = 2*num_basis_funcs * The basis must NOT be normalized. */ int *atomic_number; /**< atomic numbers (chem. element) of atoms in basis set */ - int *angular_momentum; /**< 3 ints per wave function coefficient do describe the + int *angular_momentum; /**< 3 ints per wave function coefficient do describe the * cartesian components of the angular momentum. - * E.g. S={0 0 0}, Px={1 0 0}, Dxy={1 1 0}, or Fyyz={0 2 1}. + * E.g. S={0 0 0}, Px={1 0 0}, Dxy={1 1 0}, or Fyyz={0 2 1}. */ int *shell_types; /**< type for each shell in basis */ } molfile_qm_basis_t; @@ -409,9 +409,9 @@ enum molfile_qm_wavefunc_type { MOLFILE_WAVE_MCSCFNAT, MOLFILE_WAVE_MCSCFOPT, MOLFILE_WAVE_CINATUR, MOLFILE_WAVE_PIPEK, MOLFILE_WAVE_BOYS, MOLFILE_WAVE_RUEDEN, - MOLFILE_WAVE_NAO, MOLFILE_WAVE_PNAO, MOLFILE_WAVE_NHO, - MOLFILE_WAVE_PNHO, MOLFILE_WAVE_NBO, MOLFILE_WAVE_PNBO, - MOLFILE_WAVE_PNLMO, MOLFILE_WAVE_NLMO, MOLFILE_WAVE_MOAO, + MOLFILE_WAVE_NAO, MOLFILE_WAVE_PNAO, MOLFILE_WAVE_NHO, + MOLFILE_WAVE_PNHO, MOLFILE_WAVE_NBO, MOLFILE_WAVE_PNBO, + MOLFILE_WAVE_PNLMO, MOLFILE_WAVE_NLMO, MOLFILE_WAVE_MOAO, MOLFILE_WAVE_NATO, MOLFILE_WAVE_UNKNOWN }; @@ -442,7 +442,7 @@ typedef struct molfile_qm_timestep_metadata { int has_orben_per_wavef[MOLFILE_MAXWAVEPERTS]; /**< orbital energy flags */ int has_occup_per_wavef[MOLFILE_MAXWAVEPERTS]; /**< orbital occupancy flags */ int num_wavef ; /**< # wavefunctions in this ts */ - int wavef_size; /**< size of one wavefunction + int wavef_size; /**< size of one wavefunction * (# of gaussian basis fctns) */ int num_charge_sets; /**< # of charge values per atom */ } molfile_qm_timestep_metadata_t; @@ -498,14 +498,14 @@ typedef struct { * from graphics file reader plugins. */ enum molfile_graphics_type { - MOLFILE_POINT, MOLFILE_TRIANGLE, MOLFILE_TRINORM, MOLFILE_NORMS, - MOLFILE_LINE, MOLFILE_CYLINDER, MOLFILE_CAPCYL, MOLFILE_CONE, + MOLFILE_POINT, MOLFILE_TRIANGLE, MOLFILE_TRINORM, MOLFILE_NORMS, + MOLFILE_LINE, MOLFILE_CYLINDER, MOLFILE_CAPCYL, MOLFILE_CONE, MOLFILE_SPHERE, MOLFILE_TEXT, MOLFILE_COLOR, MOLFILE_TRICOLOR }; /** * Individual graphics object/element data - */ + */ typedef struct { int type; /* One of molfile_graphics_type */ int style; /* A general style parameter */ @@ -521,13 +521,13 @@ typedef struct { type data style size ---- ---- ----- ---- point x, y, z pixel size -triangle x1,y1,z1,x2,y2,z2,x3,y3,z3 -trinorm x1,y1,z1,x2,y2,z2,x3,y3,z3 +triangle x1,y1,z1,x2,y2,z2,x3,y3,z3 +trinorm x1,y1,z1,x2,y2,z2,x3,y3,z3 the next array element must be NORMS -tricolor x1,y1,z1,x2,y2,z2,x3,y3,z3 +tricolor x1,y1,z1,x2,y2,z2,x3,y3,z3 the next array elements must be NORMS the following element must be COLOR, with three RGB triples -norms x1,y1,z1,x2,y2,z2,x3,y3,z3 +norms x1,y1,z1,x2,y2,z2,x3,y3,z3 line x1,y1,z1,x2,y2,z2 0=solid pixel width 1=stippled cylinder x1,y1,z1,x2,y2,z2 resolution radius @@ -541,41 +541,41 @@ color r, g, b /** * Main file reader API. Any function in this struct may be NULL * if not implemented by the plugin; the application checks this to determine - * what functionality is present in the plugin. - */ + * what functionality is present in the plugin. + */ typedef struct { /** - * Required header + * Required header */ vmdplugin_HEAD /** - * Filename extension for this file type. May be NULL if no filename + * Filename extension for this file type. May be NULL if no filename * extension exists and/or is known. For file types that match several * common extensions, list them in a comma separated list such as: * "pdb,ent,foo,bar,baz,ban" * The comma separated list will be expanded when filename extension matching * is performed. If multiple plugins solicit the same filename extensions, - * the one that lists the extension earliest in its list is selected. In the + * the one that lists the extension earliest in its list is selected. In the * case of a "tie", the first one tried/checked "wins". */ const char *filename_extension; /** * Try to open the file for reading. Return an opaque handle, or NULL on - * failure. Set the number of atoms; if the number of atoms cannot be - * determined, set natoms to MOLFILE_NUMATOMS_UNKNOWN. + * failure. Set the number of atoms; if the number of atoms cannot be + * determined, set natoms to MOLFILE_NUMATOMS_UNKNOWN. * Filetype should be the name under which this plugin was registered; * this is provided so that plugins can provide the same function pointer * to handle multiple file types. */ - void *(* open_file_read)(const char *filepath, const char *filetype, + void *(* open_file_read)(const char *filepath, const char *filetype, int *natoms); - + /** * Read molecular structure from the given file handle. atoms is allocated * by the caller and points to space for natoms. - * On success, place atom information in the passed-in pointer. + * On success, place atom information in the passed-in pointer. * optflags specifies which optional fields in the atoms will be set by * the plugin. */ @@ -587,15 +587,15 @@ typedef struct { * Each unique bond should be specified only once, so file formats that list * bonds twice will need post-processing before the results are returned to * the caller. - * If the plugin provides bond information, but the file loaded doesn't + * If the plugin provides bond information, but the file loaded doesn't * actually contain any bond info, the nbonds parameter should be * set to 0 and from/to should be set to NULL to indicate that no bond * information was actually present, and automatic bond search should be - * performed. + * performed. * * If the plugin provides bond order information, the bondorder array * will contain the bond order for each from/to pair. If not, the bondorder - * pointer should be set to NULL, in which case the caller will provide a + * pointer should be set to NULL, in which case the caller will provide a * default bond order value of 1.0. * * If the plugin provides bond type information, the bondtype array @@ -606,27 +606,27 @@ typedef struct { * and consistency checking. * * These arrays must be freed by the plugin in the close_file_read function. - * This function can be called only after read_structure(). - * Return MOLFILE_SUCCESS if no errors occur. + * This function can be called only after read_structure(). + * Return MOLFILE_SUCCESS if no errors occur. */ #if vmdplugin_ABIVERSION > 14 - int (*read_bonds)(void *, int *nbonds, int **from, int **to, float **bondorder, + int (*read_bonds)(void *, int *nbonds, int **from, int **to, float **bondorder, int **bondtype, int *nbondtypes, char ***bondtypename); #else int (*read_bonds)(void *, int *nbonds, int **from, int **to, float **bondorder); #endif /** - * XXX this function will be augmented and possibly superceded by a + * XXX this function will be augmented and possibly superceded by a * new QM-capable version named read_timestep(), when finished. * - * Read the next timestep from the file. Return MOLFILE_SUCCESS, or - * MOLFILE_EOF on EOF. If the molfile_timestep_t argument is NULL, then - * the frame should be skipped. Otherwise, the application must prepare - * molfile_timestep_t by allocating space in coords for the corresponding - * number of coordinates. - * The natoms parameter exists because some coordinate file formats - * (like CRD) cannot determine for themselves how many atoms are in a + * Read the next timestep from the file. Return MOLFILE_SUCCESS, or + * MOLFILE_EOF on EOF. If the molfile_timestep_t argument is NULL, then + * the frame should be skipped. Otherwise, the application must prepare + * molfile_timestep_t by allocating space in coords for the corresponding + * number of coordinates. + * The natoms parameter exists because some coordinate file formats + * (like CRD) cannot determine for themselves how many atoms are in a * timestep; the app must therefore obtain this information elsewhere * and provide it to the plugin. */ @@ -636,16 +636,16 @@ typedef struct { * Close the file and release all data. The handle cannot be reused. */ void (* close_file_read)(void *); - + /** * Open a coordinate file for writing using the given header information. * Return an opaque handle, or NULL on failure. The application must - * specify the number of atoms to be written. + * specify the number of atoms to be written. * filetype should be the name under which this plugin was registered. */ - void *(* open_file_write)(const char *filepath, const char *filetype, + void *(* open_file_write)(const char *filepath, const char *filetype, int natoms); - + /** * Write structure information. Return success. */ @@ -653,12 +653,12 @@ typedef struct { /** * Write a timestep to the coordinate file. Return MOLFILE_SUCCESS if no - * errors occur. If the file contains structure information in each - * timestep (like a multi-entry PDB), it will have to cache the information + * errors occur. If the file contains structure information in each + * timestep (like a multi-entry PDB), it will have to cache the information * from the initial calls from write_structure. */ int (* write_timestep)(void *, const molfile_timestep_t *); - + /** * Close the file and release all data. The handle cannot be reused. */ @@ -671,24 +671,24 @@ typedef struct { * the plugin and should be freed by close_file_read(). The application * may call this function any number of times. */ - int (* read_volumetric_metadata)(void *, int *nsets, + int (* read_volumetric_metadata)(void *, int *nsets, molfile_volumetric_t **metadata); - /** - * Read the specified volumetric data set into the space pointed to by - * datablock. The set is specified with a zero-based index. The space + /** + * Read the specified volumetric data set into the space pointed to by + * datablock. The set is specified with a zero-based index. The space * allocated for the datablock must be equal to - * xsize * ysize * zsize. No space will be allocated for colorblock + * xsize * ysize * zsize. No space will be allocated for colorblock * unless has_color is nonzero; in that case, colorblock should be * filled in with three RGB floats per datapoint. */ - int (* read_volumetric_data)(void *, int set, float *datablock, + int (* read_volumetric_data)(void *, int set, float *datablock, float *colorblock); /** * Read raw graphics data stored in this file. Return the number of data - * elements and the data itself as an array of molfile_graphics_t in the - * pointer provided by the application. The plugin is responsible for + * elements and the data itself as an array of molfile_graphics_t in the + * pointer provided by the application. The plugin is responsible for * freeing the data when the file is closed. */ int (* read_rawgraphics)(void *, int *nelem, const molfile_graphics_t **data); @@ -698,19 +698,19 @@ typedef struct { * came from, what the accession code for the database is, textual remarks * and other notes pertaining to the contained structure/trajectory/volume * and anything else that's informative at the whole file level. - */ + */ int (* read_molecule_metadata)(void *, molfile_metadata_t **metadata); - + /** * Write bond information for the molecule. The arrays from * and to point to the (one-based) indices of bonded atoms. - * Each unique bond will be specified only once by the caller. - * File formats that list bonds twice will need to emit both the + * Each unique bond will be specified only once by the caller. + * File formats that list bonds twice will need to emit both the * from/to and to/from versions of each. - * This function must be called before write_structure(). + * This function must be called before write_structure(). * * Like the read_bonds() routine, the bondorder pointer is set to NULL - * if the caller doesn't have such information, in which case the + * if the caller doesn't have such information, in which case the * plugin should assume a bond order of 1.0 if the file format requires * bond order information. * @@ -721,10 +721,10 @@ typedef struct { * scheme is different from the index numbers. * if the pointers are set to NULL, then this information is not available. * bondtypenames can only be used of bondtypes is also given. - * Return MOLFILE_SUCCESS if no errors occur. + * Return MOLFILE_SUCCESS if no errors occur. */ #if vmdplugin_ABIVERSION > 14 - int (* write_bonds)(void *, int nbonds, int *from, int *to, float *bondorder, + int (* write_bonds)(void *, int nbonds, int *from, int *to, float *bondorder, int *bondtype, int nbondtypes, char **bondtypename); #else int (* write_bonds)(void *, int nbonds, int *from, int *to, float *bondorder); @@ -732,9 +732,9 @@ typedef struct { #if vmdplugin_ABIVERSION > 9 /** - * Write the specified volumetric data set into the space pointed to by + * Write the specified volumetric data set into the space pointed to by * datablock. The * allocated for the datablock must be equal to - * xsize * ysize * zsize. No space will be allocated for colorblock + * xsize * ysize * zsize. No space will be allocated for colorblock * unless has_color is nonzero; in that case, colorblock should be * filled in with three RGB floats per datapoint. */ @@ -742,30 +742,30 @@ typedef struct { float *datablock, float *colorblock); #if vmdplugin_ABIVERSION > 15 - /** + /** * Read in Angles, Dihedrals, Impropers, and Cross Terms and optionally types. - * (Cross terms pertain to the CHARMM/NAMD CMAP feature) + * (Cross terms pertain to the CHARMM/NAMD CMAP feature) */ int (* read_angles)(void *handle, int *numangles, int **angles, int **angletypes, int *numangletypes, char ***angletypenames, int *numdihedrals, int **dihedrals, int **dihedraltypes, int *numdihedraltypes, - char ***dihedraltypenames, int *numimpropers, int **impropers, + char ***dihedraltypenames, int *numimpropers, int **impropers, int **impropertypes, int *numimpropertypes, char ***impropertypenames, int *numcterms, int **cterms, int *ctermcols, int *ctermrows); - /** + /** * Write out Angles, Dihedrals, Impropers, and Cross Terms - * (Cross terms pertain to the CHARMM/NAMD CMAP feature) + * (Cross terms pertain to the CHARMM/NAMD CMAP feature) */ int (* write_angles)(void *handle, int numangles, const int *angles, const int *angletypes, int numangletypes, const char **angletypenames, int numdihedrals, const int *dihedrals, const int *dihedraltypes, int numdihedraltypes, - const char **dihedraltypenames, int numimpropers, + const char **dihedraltypenames, int numimpropers, const int *impropers, const int *impropertypes, int numimpropertypes, - const char **impropertypenames, int numcterms, const int *cterms, + const char **impropertypenames, int numcterms, const int *cterms, int ctermcols, int ctermrows); #else - /** + /** * Read in Angles, Dihedrals, Impropers, and Cross Terms * Forces are in Kcal/mol * (Cross terms pertain to the CHARMM/NAMD CMAP feature, forces are given @@ -775,10 +775,10 @@ typedef struct { int *numangles, int **angles, double **angleforces, int *numdihedrals, int **dihedrals, double **dihedralforces, int *numimpropers, int **impropers, double **improperforces, - int *numcterms, int **cterms, + int *numcterms, int **cterms, int *ctermcols, int *ctermrows, double **ctermforces); - /** + /** * Write out Angles, Dihedrals, Impropers, and Cross Terms * Forces are in Kcal/mol * (Cross terms pertain to the CHARMM/NAMD CMAP feature, forces are given @@ -788,7 +788,7 @@ typedef struct { int numangles, const int *angles, const double *angleforces, int numdihedrals, const int *dihedrals, const double *dihedralforces, int numimpropers, const int *impropers, const double *improperforces, - int numcterms, const int *cterms, + int numcterms, const int *cterms, int ctermcols, int ctermrows, const double *ctermforces); #endif @@ -798,7 +798,7 @@ typedef struct { * QM datasets in this file. * * The metadata are the sizes of the QM related data structure - * arrays that will be populated by the plugin when + * arrays that will be populated by the plugin when * read_qm_rundata() is called. Since the allocation of these * arrays is done by VMD rather than the plugin, VMD needs to * know the sizes beforehand. Consequently read_qm_metadata() @@ -811,7 +811,7 @@ typedef struct { * Read timestep independent QM data. * * Typical data that are defined only once per trajectory are - * general info about the calculation (such as the used method), + * general info about the calculation (such as the used method), * the basis set and normal modes. * The data structures to be populated must have been allocated * before by VMD according to sizes obtained through @@ -821,19 +821,19 @@ typedef struct { /** - * Read the next timestep from the file. Return MOLFILE_SUCCESS, or + * Read the next timestep from the file. Return MOLFILE_SUCCESS, or * MOLFILE_EOF on EOF. If the molfile_timestep_t or molfile_qm_metadata_t - * arguments are NULL, then the coordinate or qm data should be skipped. - * Otherwise, the application must prepare molfile_timestep_t and - * molfile_qm_timestep_t by allocating space for the corresponding + * arguments are NULL, then the coordinate or qm data should be skipped. + * Otherwise, the application must prepare molfile_timestep_t and + * molfile_qm_timestep_t by allocating space for the corresponding * number of coordinates, orbital wavefunction coefficients, etc. - * Since it is common for users to want to load only the final timestep + * Since it is common for users to want to load only the final timestep * data from a QM run, the application may provide any combination of - * valid, or NULL pointers for the molfile_timestep_t and + * valid, or NULL pointers for the molfile_timestep_t and * molfile_qm_timestep_t parameters, depending on what information the * user is interested in. - * The natoms and qm metadata parameters exist because some file formats - * cannot determine for themselves how many atoms etc are in a + * The natoms and qm metadata parameters exist because some file formats + * cannot determine for themselves how many atoms etc are in a * timestep; the app must therefore obtain this information elsewhere * and provide it to the plugin. */ @@ -888,4 +888,3 @@ typedef struct { } molfile_plugin_t; #endif - diff --git a/src/USER-MOLFILE/reader_molfile.cpp b/src/USER-MOLFILE/reader_molfile.cpp index 8f4522db67..a1156efc67 100644 --- a/src/USER-MOLFILE/reader_molfile.cpp +++ b/src/USER-MOLFILE/reader_molfile.cpp @@ -36,7 +36,7 @@ enum{ID,TYPE,X,Y,Z,VX,VY,VZ}; // true if the difference between two floats is "small". // cannot use fabsf() since it is not fully portable. -static bool is_smalldiff(const float &val1, const float &val2) +static bool is_smalldiff(const float &val1, const float &val2) { return (fabs(static_cast(val1-val2)) < SMALL); } @@ -234,7 +234,7 @@ bigint ReaderMolfile::read_header(double box[3][3], int &triclinic, const double alpha = static_cast(cell[3]); const double beta = static_cast(cell[4]); const double gamma = static_cast(cell[5]); - + const double lx = la; const double xy = lb * cos(gamma/90.0*MY_PI2); const double xz = lc * cos(beta/90.0*MY_PI2); @@ -279,7 +279,7 @@ bigint ReaderMolfile::read_header(double box[3][3], int &triclinic, memory->create(fieldindex,nfield,"read_dump:fieldindex"); - // we know nothing about the scaling style of coordinates, + // we know nothing about the scaling style of coordinates, // so the caller has to set the proper flag. xflag = scaledflag; yflag = scaledflag; diff --git a/src/USER-MOLFILE/vmdplugin.h b/src/USER-MOLFILE/vmdplugin.h index f58a5e435b..37299408fe 100644 --- a/src/USER-MOLFILE/vmdplugin.h +++ b/src/USER-MOLFILE/vmdplugin.h @@ -17,20 +17,20 @@ /** @file * This header must be included by every VMD plugin library. It defines the - * API for every plugin so that VMD can organize the plugins it finds. + * API for every plugin so that VMD can organize the plugins it finds. */ #ifndef VMD_PLUGIN_H #define VMD_PLUGIN_H -/* +/* * Preprocessor tricks to make it easier for us to redefine the names of * functions when building static plugins. */ #if !defined(VMDPLUGIN) -/** - * macro defining VMDPLUGIN if it hasn't already been set to the name of +/** + * macro defining VMDPLUGIN if it hasn't already been set to the name of * a static plugin that is being compiled. This is the catch-all case. */ #define VMDPLUGIN vmdplugin @@ -38,11 +38,11 @@ /** concatenation macro, joins args x and y together as a single string */ #define xcat(x, y) cat(x, y) /** concatenation macro, joins args x and y together as a single string */ -#define cat(x, y) x ## y +#define cat(x, y) x ## y /* - * macros to correctly define plugin function names depending on whether - * the plugin is being compiled for static linkage or dynamic loading. + * macros to correctly define plugin function names depending on whether + * the plugin is being compiled for static linkage or dynamic loading. * When compiled for static linkage, each plugin needs to have unique * function names for all of its entry points. When compiled for dynamic * loading, the plugins must name their entry points consistently so that @@ -59,13 +59,13 @@ /** "WIN32" is defined on both WIN32 and WIN64 platforms... */ -#if (defined(WIN32)) +#if (defined(WIN32)) #define WIN32_LEAN_AND_MEAN #include #if !defined(STATIC_PLUGIN) #if defined(VMDPLUGIN_EXPORTS) -/** +/** * Only define DllMain for plugins, not in VMD or in statically linked plugins * VMDPLUGIN_EXPORTS is only defined when compiling dynamically loaded plugins */ @@ -86,7 +86,7 @@ BOOL APIENTRY DllMain( HANDLE hModule, #endif /* ! STATIC_PLUGIN */ #else /** If we're not compiling on Windows, then this macro is defined empty */ -#define VMDPLUGIN_API +#define VMDPLUGIN_API #endif /** define plugin linkage correctly for both C and C++ based plugins */ @@ -96,13 +96,13 @@ BOOL APIENTRY DllMain( HANDLE hModule, #define VMDPLUGIN_EXTERN extern VMDPLUGIN_API #endif /* __cplusplus */ -/* - * Plugin API functions start here +/* + * Plugin API functions start here */ -/** - * Init routine: called the first time the library is loaded by the +/** + * Init routine: called the first time the library is loaded by the * application and before any other API functions are referenced. * Return 0 on success. */ @@ -110,15 +110,15 @@ VMDPLUGIN_EXTERN int VMDPLUGIN_init(void); /** * Macro for creating a struct header used in all plugin structures. - * - * This header should be placed at the top of every plugin API definition + * + * This header should be placed at the top of every plugin API definition * so that it can be treated as a subtype of the base plugin type. * * abiversion: Defines the ABI for the base plugin type (not for other plugins) * type: A string descriptor of the plugin type. * name: A name for the plugin. * author: A string identifier, possibly including newlines. - * Major and minor version. + * Major and minor version. * is_reentrant: Whether this library can be run concurrently with itself. */ #define vmdplugin_HEAD \ @@ -129,12 +129,12 @@ VMDPLUGIN_EXTERN int VMDPLUGIN_init(void); const char *author; \ int majorv; \ int minorv; \ - int is_reentrant; + int is_reentrant; -/** +/** * Typedef for generic plugin header, individual plugins can - * make their own structures as long as the header info remains - * the same as the generic plugin header, most easily done by + * make their own structures as long as the header info remains + * the same as the generic plugin header, most easily done by * using the vmdplugin_HEAD macro. */ typedef struct { @@ -158,7 +158,7 @@ typedef struct { #define VMDPLUGIN_ERROR -1 /*@}*/ -/** +/** * Function pointer typedef for register callback functions */ typedef int (*vmdplugin_register_cb)(void *, vmdplugin_t *); @@ -175,16 +175,16 @@ typedef int (*vmdplugin_register_cb)(void *, vmdplugin_t *); VMDPLUGIN_EXTERN int VMDPLUGIN_register(void *, vmdplugin_register_cb); /** - * Allow the library to register Tcl extensions. + * Allow the library to register Tcl extensions. * This API is optional; if found by dlopen, it will be called after first - * calling init and register. + * calling init and register. */ -VMDPLUGIN_EXTERN int VMDPLUGIN_register_tcl(void *, void *tcl_interp, +VMDPLUGIN_EXTERN int VMDPLUGIN_register_tcl(void *, void *tcl_interp, vmdplugin_register_cb); /** - * The Fini method is called when the application will no longer use - * any plugins in the library. + * The Fini method is called when the application will no longer use + * any plugins in the library. */ VMDPLUGIN_EXTERN int VMDPLUGIN_fini(void); diff --git a/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp b/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp index 6ea5f6c2c1..97f7715bfb 100644 --- a/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp +++ b/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp @@ -109,7 +109,7 @@ void BondHarmonicShiftCutOMP::eval(int nfrom, int nto, ThrData * const thr) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - if (EFLAG) + if (EFLAG) ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) ); // apply force to each of 2 atoms diff --git a/src/USER-OMP/bond_harmonic_shift_omp.cpp b/src/USER-OMP/bond_harmonic_shift_omp.cpp index 2f6f1d37b9..49d985b871 100644 --- a/src/USER-OMP/bond_harmonic_shift_omp.cpp +++ b/src/USER-OMP/bond_harmonic_shift_omp.cpp @@ -106,7 +106,7 @@ void BondHarmonicShiftOMP::eval(int nfrom, int nto, ThrData * const thr) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - if (EFLAG) + if (EFLAG) ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) ); // apply force to each of 2 atoms diff --git a/src/USER-OMP/fix_qeq_comb_omp.cpp b/src/USER-OMP/fix_qeq_comb_omp.cpp index dbf073db53..18b9b5305c 100644 --- a/src/USER-OMP/fix_qeq_comb_omp.cpp +++ b/src/USER-OMP/fix_qeq_comb_omp.cpp @@ -109,7 +109,7 @@ void FixQEQCombOMP::post_force(int vflag) // more loops for first-time charge equilibrium - iloop = 0; + iloop = 0; if (firstflag) loopmax = 500; else loopmax = 200; @@ -144,7 +144,7 @@ void FixQEQCombOMP::post_force(int vflag) i = ilist[ii]; if (mask[i] & groupbit) { q1[i] += qf[i]*dtq2 - heatpq*q1[i]; - q[i] += q1[i]; + q[i] += q1[i]; } } comm->forward_comm_fix(this); @@ -152,7 +152,7 @@ void FixQEQCombOMP::post_force(int vflag) if(comb) enegtot = comb->yasu_char(qf,igroup); enegtot /= ngroup; enegchk = enegmax = 0.0; - + for (ii = 0; ii < inum ; ii++) { i = ilist[ii]; if (mask[i] & groupbit) { @@ -173,15 +173,15 @@ void FixQEQCombOMP::post_force(int vflag) if (me == 0 && fp) fprintf(fp," iteration: %d, enegtot %.6g, " "enegmax %.6g, fq deviation: %.6g\n", - iloop,enegtot,enegmax,enegchk); - + iloop,enegtot,enegmax,enegchk); + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; if (mask[i] & groupbit) - q1[i] += qf[i]*dtq2 - heatpq*q1[i]; + q1[i] += qf[i]*dtq2 - heatpq*q1[i]; } - } - + } + if (me == 0 && fp) { if (iloop == loopmax) fprintf(fp,"Charges did not converge in %d iterations\n",iloop); diff --git a/src/USER-OMP/neigh_half_bin_omp.cpp b/src/USER-OMP/neigh_half_bin_omp.cpp index 45dd79d092..ba0a180410 100644 --- a/src/USER-OMP/neigh_half_bin_omp.cpp +++ b/src/USER-OMP/neigh_half_bin_omp.cpp @@ -222,7 +222,7 @@ void Neighbor::half_bin_no_newton_ghost_omp(NeighList *list) dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; - + if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); diff --git a/src/USER-OMP/neigh_half_nsq_omp.cpp b/src/USER-OMP/neigh_half_nsq_omp.cpp index 44b777273e..6a8bee491a 100644 --- a/src/USER-OMP/neigh_half_nsq_omp.cpp +++ b/src/USER-OMP/neigh_half_nsq_omp.cpp @@ -193,12 +193,12 @@ void Neighbor::half_nsq_no_newton_ghost_omp(NeighList *list) if (includegroup && !(mask[j] & bitmask)) continue; jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; - + if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); @@ -215,16 +215,16 @@ void Neighbor::half_nsq_no_newton_ghost_omp(NeighList *list) if (includegroup && !(mask[j] & bitmask)) continue; jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; - + if (rsq <= cutneighsq[itype][jtype]) neighptr[n++] = j; } } - + ilist[i] = i; firstneigh[i] = neighptr; numneigh[i] = n; diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp index becd5ffc5d..18fba083db 100644 --- a/src/USER-REAXC/fix_qeq_reax.cpp +++ b/src/USER-REAXC/fix_qeq_reax.cpp @@ -5,7 +5,7 @@ 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 + 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. @@ -53,11 +53,11 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ -FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) : +FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - if (narg != 8) error->all(FLERR,"Illegal fix qeq/reax command"); - + if (narg != 8) error->all(FLERR,"Illegal fix qeq/reax command"); + nevery = atoi(arg[3]); swa = atof(arg[4]); swb = atof(arg[5]); @@ -73,7 +73,7 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) : s = NULL; t = NULL; nprev = 5; - + Hdia_inv = NULL; b_s = NULL; b_t = NULL; @@ -111,7 +111,7 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) : FixQEqReax::~FixQEqReax() { // unregister callbacks to this fix from Atom class - + atom->delete_callback(id,0); memory->destroy(s_hist); @@ -153,7 +153,7 @@ void FixQEqReax::pertype_parameters(char *arg) gamma = (double *) pair->extract("gamma",tmp); if (chi == NULL || eta == NULL || gamma == NULL) error->all(FLERR, - "Fix qeq/reax could not extract params from pair reax/c"); + "Fix qeq/reax could not extract params from pair reax/c"); return; } @@ -171,11 +171,11 @@ void FixQEqReax::pertype_parameters(char *arg) if (comm->me == 0) { if ((pf = fopen(arg,"r")) == NULL) error->one(FLERR,"Fix qeq/reax parameter file could not be found"); - + for (i = 1; i <= ntypes && !feof(pf); i++) { fscanf(pf,"%d %lg %lg %lg",&itype,&v1,&v2,&v3); if (itype < 1 || itype > ntypes) - error->one(FLERR,"Fix qeq/reax invalid atom type in param file"); + error->one(FLERR,"Fix qeq/reax invalid atom type in param file"); chi[itype] = v1; eta[itype] = v2; gamma[itype] = v3; @@ -255,7 +255,7 @@ void FixQEqReax::allocate_matrix() m += list->numneigh[i]; } m_cap = MAX( (int)(m * SAFE_ZONE), MIN_CAP * MIN_NBRS ); - + H.n = n_cap; H.m = m_cap; memory->create(H.firstnbr,n_cap,"qeq:H.firstnbr"); @@ -271,7 +271,7 @@ void FixQEqReax::deallocate_matrix() memory->destroy( H.firstnbr ); memory->destroy( H.numnbrs ); memory->destroy( H.jlist ); - memory->destroy( H.val ); + memory->destroy( H.val ); } /* ---------------------------------------------------------------------- */ @@ -287,7 +287,7 @@ void FixQEqReax::reallocate_matrix() void FixQEqReax::init() { if (!atom->q_flag) error->all(FLERR,"Fix qeq/reax requires atom attribute q"); - + // need a half neighbor list w/ Newton off and ghost neighbors // built whenever re-neighboring occurs @@ -320,7 +320,7 @@ void FixQEqReax::init_shielding() ntypes = atom->ntypes; memory->create(shld,ntypes+1,ntypes+1,"qeq:shileding"); - + for( i = 1; i <= ntypes; ++i ) for( j = 1; j <= ntypes; ++j ) shld[i][j] = pow( gamma[i] * gamma[j], -1.5 ); @@ -353,7 +353,7 @@ void FixQEqReax::init_taper() Tap[2] =-210.0 * (swa3*swb2 + swa2*swb3) / d7; Tap[1] = 140.0 * swa3 * swb3 / d7; Tap[0] = (-35.0*swa3*swb2*swb2 + 21.0*swa2*swb3*swb2 + - 7.0*swa*swb3*swb3 + swb3*swb3*swb ) / d7; + 7.0*swa*swb3*swb3 + swb3*swb3*swb ) / d7; } /* ---------------------------------------------------------------------- */ @@ -366,7 +366,7 @@ void FixQEqReax::setup_pre_force(int vflag) allocate_matrix(); pre_force(vflag); -} +} /* ---------------------------------------------------------------------- */ @@ -374,14 +374,14 @@ void FixQEqReax::setup_pre_force_respa(int vflag, int ilevel) { if (ilevel < nlevels_respa-1) return; setup_pre_force(vflag); -} +} /* ---------------------------------------------------------------------- */ void FixQEqReax::min_setup_pre_force(int vflag) { setup_pre_force(vflag); -} +} /* ---------------------------------------------------------------------- */ @@ -389,7 +389,7 @@ void FixQEqReax::init_storage() { N = atom->nlocal + atom->nghost; for( int i = 0; i < N; i++ ) { - Hdia_inv[i] = 1. / eta[atom->type[i]]; + Hdia_inv[i] = 1. / eta[atom->type[i]]; b_s[i] = -chi[atom->type[i]]; b_t[i] = -1.0; b_prc[i] = 0; @@ -414,9 +414,9 @@ void FixQEqReax::pre_force(int vflag) if( atom->nmax > nmax ) reallocate_storage(); if( n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE ) reallocate_matrix(); - + init_matvec(); - matvecs = CG(b_s, s); // CG on s - parallel + matvecs = CG(b_s, s); // CG on s - parallel matvecs += CG(b_t, t); // CG on t - parallel calculate_Q(); @@ -458,7 +458,7 @@ void FixQEqReax::init_matvec() //t[i] = 2 * t_hist[i][0] - t_hist[i][1]; /* quadratic extrapolation for s & t from previous solutions */ - //s[i] = s_hist[i][2] + 3 * ( s_hist[i][0] - s_hist[i][1] ); + //s[i] = s_hist[i][2] + 3 * ( s_hist[i][0] - s_hist[i][1] ); t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] ); /* cubic extrapolation for s & t from previous solutions */ @@ -490,7 +490,7 @@ void FixQEqReax::compute_H() ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - + // fill in the H matrix m_fill = 0; r_sqr = 0; @@ -499,20 +499,20 @@ void FixQEqReax::compute_H() jlist = firstneigh[i]; jnum = numneigh[i]; H.firstnbr[i] = m_fill; - + for( jj = 0; jj < jnum; jj++ ) { j = jlist[jj]; - + dx = x[j][0] - x[i][0]; dy = x[j][1] - x[i][1]; dz = x[j][2] - x[i][2]; r_sqr = SQR(dx) + SQR(dy) + SQR(dz); - + flag = 0; if (r_sqr <= SQR(swb)) { if (j < n) flag = 1; else if (tag[i] < tag[j]) flag = 1; - else if (tag[i] == tag[j]) { + else if (tag[i] == tag[j]) { if (dz > SMALL) flag = 1; else if (fabs(dz) < SMALL) { if (dy > SMALL) flag = 1; @@ -521,21 +521,21 @@ void FixQEqReax::compute_H() } } } - + if( flag ) { - H.jlist[m_fill] = j; - H.val[m_fill] = calculate_H( sqrt(r_sqr), shld[type[i]][type[j]] ); - m_fill++; + H.jlist[m_fill] = j; + H.val[m_fill] = calculate_H( sqrt(r_sqr), shld[type[i]][type[j]] ); + m_fill++; } } - + H.numnbrs[i] = m_fill - H.firstnbr[i]; } if (m_fill >= H.m) { char str[128]; sprintf(str,"H matrix size has been exceeded: m_fill=%d H.m=%d\n", - m_fill, H.m ); + m_fill, H.m ); error->warning(FLERR,str); error->all(FLERR,"Fix qeq/reax has insufficient QEq matrix size"); } @@ -585,21 +585,21 @@ int FixQEqReax::CG( double *b, double *x ) comm->forward_comm_fix(this); //Dist_vector( d ); sparse_matvec( &H, d, q ); comm->reverse_comm_fix(this); //Coll_vector( q ); - + tmp = parallel_dot( d, q, n ); alpha = sig_new / tmp; // comm->me, i, parallel_norm( d, n ), parallel_norm( q, n ), tmp ); - + vector_add( x, alpha, d, n ); vector_add( r, -alpha, q, n ); - + // pre-conditioning for( j = 0; j < n; ++j ) p[j] = r[j] * Hdia_inv[j]; - + sig_old = sig_new; sig_new = parallel_dot( r, p, n ); - + beta = sig_new / sig_old; vector_sum( d, 1., p, beta, d, n ); @@ -621,7 +621,7 @@ void FixQEqReax::sparse_matvec( sparse_matrix *A, double *x, double *b ) b[i] = eta[ atom->type[i] ] * x[i]; for( i = n; i < N; ++i ) b[i] = 0; - + for( i = 0; i < n; ++i ) { for( itr_j=A->firstnbr[i]; itr_jfirstnbr[i]+A->numnbrs[i]; itr_j++) { j = A->jlist[itr_j]; @@ -642,10 +642,10 @@ void FixQEqReax::calculate_Q() s_sum = parallel_vector_acc( s, n ); t_sum = parallel_vector_acc( t, n); u = s_sum / t_sum; - + for( i = 0; i < n; ++i ) { q[i] = s[i] - u * t[i]; - + /* backup s & t */ for( k = 4; k > 0; --k ) { s_hist[i][k] = s_hist[i][k-1]; @@ -661,12 +661,12 @@ void FixQEqReax::calculate_Q() /* ---------------------------------------------------------------------- */ -int FixQEqReax::pack_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) +int FixQEqReax::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) { int m; - if( pack_flag == 1) + if( pack_flag == 1) for(m = 0; m < n; m++) buf[m] = d[list[m]]; else if( pack_flag == 2 ) for(m = 0; m < n; m++) buf[m] = s[list[m]]; @@ -677,23 +677,23 @@ int FixQEqReax::pack_comm(int n, int *list, double *buf, return 1; } - + /* ---------------------------------------------------------------------- */ void FixQEqReax::unpack_comm(int n, int first, double *buf) { int i, m; - - if( pack_flag == 1) + + if( pack_flag == 1) for(m = 0, i = first; m < n; m++, i++) d[i] = buf[m]; - else if( pack_flag == 2) + else if( pack_flag == 2) for(m = 0, i = first; m < n; m++, i++) s[i] = buf[m]; - else if( pack_flag == 3) + else if( pack_flag == 3) for(m = 0, i = first; m < n; m++, i++) t[i] = buf[m]; - else if( pack_flag == 4) + else if( pack_flag == 4) for(m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m]; } - + /* ---------------------------------------------------------------------- */ int FixQEqReax::pack_reverse_comm(int n, int first, double *buf) @@ -702,7 +702,7 @@ int FixQEqReax::pack_reverse_comm(int n, int first, double *buf) for(m = 0, i = first; m < n; m++, i++) buf[m] = q[i]; return 1; } - + /* ---------------------------------------------------------------------- */ void FixQEqReax::unpack_reverse_comm(int n, int *list, double *buf) @@ -722,7 +722,7 @@ double FixQEqReax::memory_usage() bytes += atom->nmax*11 * sizeof(double); // storage bytes += n_cap*2 * sizeof(int); // matrix... bytes += m_cap * sizeof(int); - bytes += m_cap * sizeof(double); + bytes += m_cap * sizeof(double); return bytes; } @@ -825,7 +825,7 @@ double FixQEqReax::parallel_vector_acc( double *v, int n ) double FixQEqReax::norm( double* v1, int k ) { double ret = 0; - + for( --k; k>=0; --k ) ret += ( v1[k] * v1[k] ); @@ -834,8 +834,8 @@ double FixQEqReax::norm( double* v1, int k ) /* ---------------------------------------------------------------------- */ -void FixQEqReax::vector_sum( double* dest, double c, double* v, - double d, double* y, int k ) +void FixQEqReax::vector_sum( double* dest, double c, double* v, + double d, double* y, int k ) { for( --k; k>=0; --k ) dest[k] = c * v[k] + d * y[k]; @@ -854,7 +854,7 @@ void FixQEqReax::vector_scale( double* dest, double c, double* v, int k ) double FixQEqReax::dot( double* v1, double* v2, int k ) { double ret = 0; - + for( --k; k>=0; --k ) ret += v1[k] * v2[k]; diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp index cddf6d8a9a..47821b1056 100644 --- a/src/USER-REAXC/pair_reax_c.cpp +++ b/src/USER-REAXC/pair_reax_c.cpp @@ -5,7 +5,7 @@ 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 + 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. @@ -62,13 +62,13 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp) system = (reax_system *) memory->smalloc(sizeof(reax_system),"reax:system"); - control = (control_params *) + control = (control_params *) memory->smalloc(sizeof(control_params),"reax:control"); data = (simulation_data *) memory->smalloc(sizeof(simulation_data),"reax:data"); workspace = (storage *) memory->smalloc(sizeof(storage),"reax:storage"); - lists = (reax_list *) + lists = (reax_list *) memory->smalloc(LIST_N * sizeof(reax_list),"reax:lists"); out_control = (output_controls *) memory->smalloc(sizeof(output_controls),"reax:out_control"); @@ -128,7 +128,7 @@ PairReaxC::~PairReaxC() DeAllocate_System( system ); } //fprintf( stderr, "4\n" ); - + memory->destroy( system ); memory->destroy( control ); memory->destroy( data ); @@ -192,7 +192,7 @@ void PairReaxC::settings(int narg, char **arg) control->hbond_cut = 7.50; control->thb_cut = 0.001; control->thb_cutsq = 0.00001; - + out_control->write_steps = 0; out_control->traj_method = 0; strcpy( out_control->traj_title, "default_title" ); @@ -248,7 +248,7 @@ void PairReaxC::coeff( int nargs, char **args ) // read ffield file Read_Force_Field(args[2], &(system->reax_param), control); - + // read args that map atom types to elements in potential file // map[i] = which element the Ith atom type is, -1 if NULL @@ -278,7 +278,7 @@ void PairReaxC::coeff( int nargs, char **args ) setflag[i][j] = 1; count++; } - + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -293,7 +293,7 @@ void PairReaxC::init_style( ) int iqeq; for (iqeq = 0; iqeq < modify->nfix; iqeq++) if (strcmp(modify->fix[iqeq]->style,"qeq/reax") == 0) break; - if (iqeq == modify->nfix && qeqflag == 1) + if (iqeq == modify->nfix && qeqflag == 1) error->all(FLERR,"Pair reax/c requires use of fix qeq/reax"); system->n = atom->nlocal; // my atoms @@ -302,8 +302,8 @@ void PairReaxC::init_style( ) system->wsize = comm->nprocs; system->big_box.V = 0; - system->big_box.box_norms[0] = 0; - system->big_box.box_norms[1] = 0; + system->big_box.box_norms[0] = 0; + system->big_box.box_norms[1] = 0; system->big_box.box_norms[2] = 0; if (atom->tag_enable == 0) @@ -348,7 +348,7 @@ void PairReaxC::setup( ) if (setup_flag == 0) { setup_flag = 1; - + int *num_bonds = fix_reax->num_bonds; int *num_hbonds = fix_reax->num_hbonds; @@ -363,15 +363,15 @@ void PairReaxC::setup( ) PreAllocate_Space( system, control, workspace, world ); write_reax_atoms(); - + int num_nbrs = estimate_reax_lists(); - if(!Make_List(system->total_cap, num_nbrs, TYP_FAR_NEIGHBOR, - lists+FAR_NBRS, world)) + if(!Make_List(system->total_cap, num_nbrs, TYP_FAR_NEIGHBOR, + lists+FAR_NBRS, world)) error->all(FLERR,"Pair reax/c problem in far neighbor list"); - + write_reax_lists(); - Initialize( system, control, data, workspace, &lists, out_control, - mpi_data, world ); + Initialize( system, control, data, workspace, &lists, out_control, + mpi_data, world ); for( int k = 0; k < system->N; ++k ) { num_bonds[k] = system->my_atoms[k].num_bonds; num_hbonds[k] = system->my_atoms[k].num_hbonds; @@ -387,7 +387,7 @@ void PairReaxC::setup( ) for(int k = oldN; k < system->N; ++k) Set_End_Index( k, Start_Index( k, lists+BONDS ), lists+BONDS ); - + // check if I need to shrink/extend my data-structs ReAllocate( system, control, data, workspace, &lists, mpi_data ); @@ -422,9 +422,9 @@ void PairReaxC::compute(int eflag, int vflag) /* if ((eflag_atom || vflag_atom) && firstwarn) { firstwarn = 0; - if (comm->me == 0) + if (comm->me == 0) error->warning(FLERR,"Pair reax/c cannot yet compute " - "per-atom energy or stress"); + "per-atom energy or stress"); } */ if (vflag_global) control->virial = 1; @@ -435,15 +435,15 @@ void PairReaxC::compute(int eflag, int vflag) system->bigN = static_cast (atom->natoms); // all atoms in the system system->big_box.V = 0; - system->big_box.box_norms[0] = 0; - system->big_box.box_norms[1] = 0; + system->big_box.box_norms[0] = 0; + system->big_box.box_norms[1] = 0; system->big_box.box_norms[2] = 0; if( comm->me == 0 ) t_start = MPI_Wtime(); // setup data structures setup(); - + Reset( system, control, data, workspace, &lists, world ); workspace->realloc.num_far = write_reax_lists(); // timing for filling in the reax lists @@ -486,7 +486,7 @@ void PairReaxC::compute(int eflag, int vflag) // Store the different parts of the energy // in a list for output by compute pair command - pvector[0] = data->my_en.e_bond; + pvector[0] = data->my_en.e_bond; pvector[1] = data->my_en.e_ov + data->my_en.e_un; pvector[2] = data->my_en.e_lp; pvector[3] = 0.0; @@ -517,11 +517,11 @@ void PairReaxC::compute(int eflag, int vflag) Output_Results( system, control, data, &lists, out_control, mpi_data ); - if(fixbond_flag) - fixbond( system, control, data, &lists, out_control, mpi_data ); + if(fixbond_flag) + fixbond( system, control, data, &lists, out_control, mpi_data ); - if(fixspecies_flag) - fixspecies( system, control, data, &lists, out_control, mpi_data ); + if(fixspecies_flag) + fixspecies( system, control, data, &lists, out_control, mpi_data ); } @@ -531,7 +531,7 @@ void PairReaxC::write_reax_atoms() { int *num_bonds = fix_reax->num_bonds; int *num_hbonds = fix_reax->num_hbonds; - + for( int i = 0; i < system->N; ++i ){ system->my_atoms[i].orig_id = atom->tag[i]; system->my_atoms[i].type = map[atom->type[i]]; @@ -556,8 +556,8 @@ void PairReaxC::get_distance( rvec xj, rvec xi, double *d_sqr, rvec *dvec ) /* ---------------------------------------------------------------------- */ -void PairReaxC::set_far_nbr( far_neighbor_data *fdest, - int j, double d, rvec dvec ) +void PairReaxC::set_far_nbr( far_neighbor_data *fdest, + int j, double d, rvec dvec ) { fdest->nbr = j; fdest->d = d; @@ -607,9 +607,9 @@ int PairReaxC::estimate_reax_lists() j = jlist[itr_j]; j &= NEIGHMASK; get_distance( x[j], x[i], &d_sqr, &dvec ); - + if( d_sqr <= SQR(control->nonb_cut) ) - ++num_nbrs; + ++num_nbrs; } } @@ -643,7 +643,7 @@ int PairReaxC::write_reax_lists() far_nbrs = lists + FAR_NBRS; far_list = far_nbrs->select.far_nbr_list; - num_nbrs = 0; + num_nbrs = 0; marked = (int*) calloc( system->N, sizeof(int) ); dist = (double*) calloc( system->N, sizeof(double) ); @@ -664,8 +664,8 @@ int PairReaxC::write_reax_lists() if( d_sqr <= (control->nonb_cut*control->nonb_cut) ){ dist[j] = sqrt( d_sqr ); - set_far_nbr( &far_list[num_nbrs], j, dist[j], dvec ); - ++num_nbrs; + set_far_nbr( &far_list[num_nbrs], j, dist[j], dvec ); + ++num_nbrs; } } Set_End_Index( i, num_nbrs, far_nbrs ); @@ -673,12 +673,12 @@ int PairReaxC::write_reax_lists() free( marked ); free( dist ); - + return num_nbrs; } - + /* ---------------------------------------------------------------------- */ - + void PairReaxC::read_reax_forces() { for( int i = 0; i < system->N; ++i ) { diff --git a/src/USER-REAXC/reaxc_nonbonded.cpp b/src/USER-REAXC/reaxc_nonbonded.cpp index 691900554c..b244954881 100644 --- a/src/USER-REAXC/reaxc_nonbonded.cpp +++ b/src/USER-REAXC/reaxc_nonbonded.cpp @@ -327,12 +327,12 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control, //fprintf(stderr, "r: %f, i: %d, base: %f, dif: %f\n", r, i, base, dif); e_vdW = ((t->vdW[r].d*dif + t->vdW[r].c)*dif + t->vdW[r].b)*dif + - t->vdW[r].a; + t->vdW[r].a; e_ele = ((t->ele[r].d*dif + t->ele[r].c)*dif + t->ele[r].b)*dif + - t->ele[r].a; + t->ele[r].a; e_ele *= system->my_atoms[i].q * system->my_atoms[j].q; - + data->my_en.e_vdW += e_vdW; data->my_en.e_ele += e_ele; diff --git a/src/USER-SPH/atom_vec_meso.cpp b/src/USER-SPH/atom_vec_meso.cpp index 6ca52f0f84..73cf49619f 100644 --- a/src/USER-SPH/atom_vec_meso.cpp +++ b/src/USER-SPH/atom_vec_meso.cpp @@ -31,7 +31,7 @@ AtomVecMeso::AtomVecMeso(LAMMPS *lmp, int narg, char **arg) : AtomVec(lmp, narg, arg) { molecular = 0; mass_type = 1; - + comm_x_only = 0; // we communicate not only x forward but also vest ... comm_f_only = 0; // we also communicate de and drho in reverse direction size_forward = 8; // 3 + rho + e + vest[3], that means we may only communicate 5 in hybrid @@ -41,7 +41,7 @@ AtomVecMeso::AtomVecMeso(LAMMPS *lmp, int narg, char **arg) : size_data_atom = 8; size_data_vel = 4; xcol_data = 6; - + atom->e_flag = 1; atom->rho_flag = 1; atom->cv_flag = 1; @@ -62,7 +62,7 @@ void AtomVecMeso::grow(int n) { atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) error->one(FLERR,"Per-processor system is too big"); - + tag = memory->grow(atom->tag, nmax, "atom:tag"); type = memory->grow(atom->type, nmax, "atom:type"); mask = memory->grow(atom->mask, nmax, "atom:mask"); @@ -70,14 +70,14 @@ void AtomVecMeso::grow(int n) { x = memory->grow(atom->x, nmax, 3, "atom:x"); v = memory->grow(atom->v, nmax, 3, "atom:v"); f = memory->grow(atom->f, nmax*comm->nthreads, 3, "atom:f"); - + rho = memory->grow(atom->rho, nmax, "atom:rho"); drho = memory->grow(atom->drho, nmax*comm->nthreads, "atom:drho"); e = memory->grow(atom->e, nmax, "atom:e"); de = memory->grow(atom->de, nmax*comm->nthreads, "atom:de"); vest = memory->grow(atom->vest, nmax, 3, "atom:vest"); cv = memory->grow(atom->cv, nmax, "atom:cv"); - + if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); @@ -117,7 +117,7 @@ void AtomVecMeso::copy(int i, int j, int delflag) { v[j][0] = v[i][0]; v[j][1] = v[i][1]; v[j][2] = v[i][2]; - + rho[j] = rho[i]; drho[j] = drho[i]; e[j] = e[i]; @@ -126,7 +126,7 @@ void AtomVecMeso::copy(int i, int j, int delflag) { vest[j][0] = vest[i][0]; vest[j][1] = vest[i][1]; vest[j][2] = vest[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); @@ -137,7 +137,7 @@ void AtomVecMeso::copy(int i, int j, int delflag) { int AtomVecMeso::pack_comm_hybrid(int n, int *list, double *buf) { //printf("in AtomVecMeso::pack_comm_hybrid\n"); int i, j, m; - + m = 0; for (i = 0; i < n; i++) { j = list[i]; @@ -155,7 +155,7 @@ int AtomVecMeso::pack_comm_hybrid(int n, int *list, double *buf) { int AtomVecMeso::unpack_comm_hybrid(int n, int first, double *buf) { //printf("in AtomVecMeso::unpack_comm_hybrid\n"); int i, m, last; - + m = 0; last = first + n; for (i = first; i < last; i++) { @@ -173,7 +173,7 @@ int AtomVecMeso::unpack_comm_hybrid(int n, int first, double *buf) { int AtomVecMeso::pack_border_hybrid(int n, int *list, double *buf) { //printf("in AtomVecMeso::pack_border_hybrid\n"); int i, j, m; - + m = 0; for (i = 0; i < n; i++) { j = list[i]; @@ -191,7 +191,7 @@ int AtomVecMeso::pack_border_hybrid(int n, int *list, double *buf) { int AtomVecMeso::unpack_border_hybrid(int n, int first, double *buf) { //printf("in AtomVecMeso::unpack_border_hybrid\n"); int i, m, last; - + m = 0; last = first + n; for (i = first; i < last; i++) { @@ -209,7 +209,7 @@ int AtomVecMeso::unpack_border_hybrid(int n, int first, double *buf) { int AtomVecMeso::pack_reverse_hybrid(int n, int first, double *buf) { //printf("in AtomVecMeso::pack_reverse_hybrid\n"); int i, m, last; - + m = 0; last = first + n; for (i = first; i < last; i++) { @@ -224,7 +224,7 @@ int AtomVecMeso::pack_reverse_hybrid(int n, int first, double *buf) { int AtomVecMeso::unpack_reverse_hybrid(int n, int *list, double *buf) { //printf("in AtomVecMeso::unpack_reverse_hybrid\n"); int i, j, m; - + m = 0; for (i = 0; i < n; i++) { j = list[i]; @@ -241,7 +241,7 @@ int AtomVecMeso::pack_comm(int n, int *list, double *buf, int pbc_flag, //printf("in AtomVecMeso::pack_comm\n"); int i, j, m; double dx, dy, dz; - + m = 0; if (pbc_flag == 0) { for (i = 0; i < n; i++) { @@ -287,7 +287,7 @@ int AtomVecMeso::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, //printf("in AtomVecMeso::pack_comm_vel\n"); int i, j, m; double dx, dy, dz; - + m = 0; if (pbc_flag == 0) { for (i = 0; i < n; i++) { @@ -337,7 +337,7 @@ int AtomVecMeso::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, void AtomVecMeso::unpack_comm(int n, int first, double *buf) { //printf("in AtomVecMeso::unpack_comm\n"); int i, m, last; - + m = 0; last = first + n; for (i = first; i < last; i++) { @@ -357,7 +357,7 @@ void AtomVecMeso::unpack_comm(int n, int first, double *buf) { void AtomVecMeso::unpack_comm_vel(int n, int first, double *buf) { //printf("in AtomVecMeso::unpack_comm_vel\n"); int i, m, last; - + m = 0; last = first + n; for (i = first; i < last; i++) { @@ -380,7 +380,7 @@ void AtomVecMeso::unpack_comm_vel(int n, int first, double *buf) { int AtomVecMeso::pack_reverse(int n, int first, double *buf) { //printf("in AtomVecMeso::pack_reverse\n"); int i, m, last; - + m = 0; last = first + n; for (i = first; i < last; i++) { @@ -398,7 +398,7 @@ int AtomVecMeso::pack_reverse(int n, int first, double *buf) { void AtomVecMeso::unpack_reverse(int n, int *list, double *buf) { //printf("in AtomVecMeso::unpack_reverse\n"); int i, j, m; - + m = 0; for (i = 0; i < n; i++) { j = list[i]; @@ -417,7 +417,7 @@ int AtomVecMeso::pack_border(int n, int *list, double *buf, int pbc_flag, //printf("in AtomVecMeso::pack_border\n"); int i, j, m; double dx, dy, dz; - + m = 0; if (pbc_flag == 0) { for (i = 0; i < n; i++) { @@ -471,7 +471,7 @@ int AtomVecMeso::pack_border_vel(int n, int *list, double *buf, int pbc_flag, //printf("in AtomVecMeso::pack_border_vel\n"); int i, j, m; double dx, dy, dz; - + m = 0; if (pbc_flag == 0) { for (i = 0; i < n; i++) { @@ -529,7 +529,7 @@ int AtomVecMeso::pack_border_vel(int n, int *list, double *buf, int pbc_flag, void AtomVecMeso::unpack_border(int n, int first, double *buf) { //printf("in AtomVecMeso::unpack_border\n"); int i, m, last; - + m = 0; last = first + n; for (i = first; i < last; i++) { @@ -555,7 +555,7 @@ void AtomVecMeso::unpack_border(int n, int first, double *buf) { void AtomVecMeso::unpack_border_vel(int n, int first, double *buf) { //printf("in AtomVecMeso::unpack_border_vel\n"); int i, m, last; - + m = 0; last = first + n; for (i = first; i < last; i++) { @@ -603,11 +603,11 @@ int AtomVecMeso::pack_exchange(int i, double *buf) { buf[m++] = vest[i][0]; buf[m++] = vest[i][1]; buf[m++] = vest[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; } @@ -619,7 +619,7 @@ int AtomVecMeso::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++]; @@ -637,12 +637,12 @@ int AtomVecMeso::unpack_exchange(double *buf) { vest[nlocal][0] = buf[m++]; vest[nlocal][1] = buf[m++]; vest[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; } @@ -654,15 +654,15 @@ int AtomVecMeso::unpack_exchange(double *buf) { int AtomVecMeso::size_restart() { int i; - + int nlocal = atom->nlocal; int n = 17 * nlocal; // 11 + rho + e + cv + vest[3] - + 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; } @@ -690,11 +690,11 @@ int AtomVecMeso::pack_restart(int i, double *buf) { buf[m++] = vest[i][0]; buf[m++] = vest[i][1]; buf[m++] = vest[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; } @@ -710,7 +710,7 @@ int AtomVecMeso::unpack_restart(double *buf) { if (atom->nextra_store) memory->grow(atom->extra, nmax, atom->nextra_store, "atom:extra"); } - + int m = 1; x[nlocal][0] = buf[m++]; x[nlocal][1] = buf[m++]; @@ -728,14 +728,14 @@ int AtomVecMeso::unpack_restart(double *buf) { vest[nlocal][0] = buf[m++]; vest[nlocal][1] = buf[m++]; vest[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; } @@ -749,14 +749,14 @@ void AtomVecMeso::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] = ((tagint) IMGMAX << IMG2BITS) | + image[nlocal] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMAX; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; @@ -769,7 +769,7 @@ void AtomVecMeso::create_atom(int itype, double *coord) { vest[nlocal][2] = 0.0; de[nlocal] = 0.0; drho[nlocal] = 0.0; - + atom->nlocal++; } @@ -782,39 +782,39 @@ void AtomVecMeso::data_atom(double *coord, tagint imagetmp, char **values) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); - + tag[nlocal] = atoi(values[0]); if (tag[nlocal] <= 0) error->one(FLERR,"Invalid atom ID in Atoms section of data file"); - + type[nlocal] = atoi(values[1]); if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) error->one(FLERR,"Invalid atom type in Atoms section of data file"); - + rho[nlocal] = atof(values[2]); e[nlocal] = atof(values[3]); cv[nlocal] = atof(values[4]); - + x[nlocal][0] = coord[0]; x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; - + //printf("rho=%f, e=%f, cv=%f, x=%f\n", rho[nlocal], e[nlocal], cv[nlocal], x[nlocal][0]); - + image[nlocal] = imagetmp; - + mask[nlocal] = 1; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; v[nlocal][2] = 0.0; - + vest[nlocal][0] = 0.0; vest[nlocal][1] = 0.0; vest[nlocal][2] = 0.0; - + de[nlocal] = 0.0; drho[nlocal] = 0.0; - + atom->nlocal++; } @@ -824,11 +824,11 @@ void AtomVecMeso::data_atom(double *coord, tagint imagetmp, char **values) { ------------------------------------------------------------------------- */ int AtomVecMeso::data_atom_hybrid(int nlocal, char **values) { - + rho[nlocal] = atof(values[0]); e[nlocal] = atof(values[1]); cv[nlocal] = atof(values[2]); - + return 3; } @@ -838,7 +838,7 @@ int AtomVecMeso::data_atom_hybrid(int nlocal, char **values) { bigint AtomVecMeso::memory_usage() { bigint bytes = 0; - + if (atom->memcheck("tag")) bytes += memory->usage(tag, nmax); if (atom->memcheck("type")) @@ -865,6 +865,6 @@ bigint AtomVecMeso::memory_usage() { bytes += memory->usage(cv, nmax); if (atom->memcheck("vest")) bytes += memory->usage(vest, nmax); - + return bytes; }