forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7067 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
90a5d2b322
commit
855e2aa353
|
|
@ -1,166 +0,0 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "lmptype.h"
|
||||
#include "mpi.h"
|
||||
#include <math.h>
|
||||
#include "fix_qeq_comb_omp.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "group.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "respa.h"
|
||||
#include "update.h"
|
||||
#include "pair_comb_omp.h"
|
||||
|
||||
#include <string.h>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixQEQCombOMP::FixQEQCombOMP(LAMMPS *lmp, int narg, char **arg)
|
||||
: FixQEQComb(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 5) error->all(FLERR,"Illegal fix qeq/comb/omp command");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEQCombOMP::init()
|
||||
{
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,"Fix qeq/comb/omp requires atom attribute q");
|
||||
|
||||
comb = (PairComb *) force->pair_match("comb/omp",1);
|
||||
if (comb == NULL)
|
||||
comb = (PairComb *) force->pair_match("comb",1);
|
||||
if (comb == NULL) error->all(FLERR,"Must use pair_style comb or comb/omp with fix qeq/comb");
|
||||
|
||||
if (strstr(update->integrate_style,"respa"))
|
||||
nlevels_respa = ((Respa *) update->integrate)->nlevels;
|
||||
|
||||
ngroup = group->count(igroup);
|
||||
if (ngroup == 0) error->all(FLERR,"Fix qeq/comb group has no atoms");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEQCombOMP::post_force(int vflag)
|
||||
{
|
||||
int i,iloop,loopmax;
|
||||
double heatpq,qmass,dtq,dtq2;
|
||||
double enegchkall,enegmaxall;
|
||||
|
||||
if (update->ntimestep % nevery) return;
|
||||
|
||||
// reallocate work arrays if necessary
|
||||
// qf = charge force
|
||||
// q1 = charge displacement
|
||||
// q2 = tmp storage of charge force for next iteration
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(qf);
|
||||
memory->destroy(q1);
|
||||
memory->destroy(q2);
|
||||
nmax = atom->nmax;
|
||||
memory->create(qf,nmax,"qeq:qf");
|
||||
memory->create(q1,nmax,"qeq:q1");
|
||||
memory->create(q2,nmax,"qeq:q2");
|
||||
vector_atom = qf;
|
||||
}
|
||||
|
||||
// more loops for first-time charge equilibrium
|
||||
|
||||
iloop = 0;
|
||||
if (firstflag) loopmax = 5000;
|
||||
else loopmax = 2000;
|
||||
|
||||
// charge-equilibration loop
|
||||
|
||||
if (me == 0 && fp)
|
||||
fprintf(fp,"Charge equilibration on step " BIGINT_FORMAT "\n",
|
||||
update->ntimestep);
|
||||
|
||||
heatpq = 0.05;
|
||||
qmass = 0.000548580;
|
||||
dtq = 0.0006;
|
||||
dtq2 = 0.5*dtq*dtq/qmass;
|
||||
|
||||
double enegchk = 0.0;
|
||||
double enegtot = 0.0;
|
||||
double enegmax = 0.0;
|
||||
|
||||
double *q = atom->q;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
q1[i] = q2[i] = qf[i] = 0.0;
|
||||
|
||||
for (iloop = 0; iloop < loopmax; iloop ++ ) {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
q1[i] += qf[i]*dtq2 - heatpq*q1[i];
|
||||
q[i] += q1[i];
|
||||
}
|
||||
|
||||
enegtot = comb->yasu_char(qf,igroup);
|
||||
enegtot /= ngroup;
|
||||
enegchk = enegmax = 0.0;
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel for private(i) default(shared)
|
||||
#endif
|
||||
for (i = 0; i < nlocal ; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
q2[i] = enegtot-qf[i];
|
||||
enegmax = MAX(enegmax,fabs(q2[i]));
|
||||
enegchk += fabs(q2[i]);
|
||||
qf[i] = q2[i];
|
||||
}
|
||||
|
||||
MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
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);
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel for private(i) default(shared)
|
||||
#endif
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -1,32 +0,0 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(qeq/comb/omp,FixQEQCombOMP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_QEQ_COMB_OMP_H
|
||||
#define LMP_FIX_QEQ_COMB_OMP_H
|
||||
|
||||
#include "fix_qeq_comb.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixQEQCombOMP : public FixQEQComb {
|
||||
public:
|
||||
FixQEQCombOMP(class LAMMPS *, int, char **);
|
||||
virtual void init();
|
||||
virtual void post_force(int);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -1,540 +0,0 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
This software is distributed under the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "pair_comb_omp.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "group.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCombOMP::PairCombOMP(LAMMPS *lmp) :
|
||||
PairComb(lmp), ThrOMP(lmp, PAIR)
|
||||
{
|
||||
respa_enable = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairCombOMP::compute(int eflag, int vflag)
|
||||
{
|
||||
if (eflag || vflag) {
|
||||
ev_setup(eflag,vflag);
|
||||
ev_setup_thr(this);
|
||||
} else evflag = vflag_fdotr = vflag_atom = 0;
|
||||
|
||||
const int nall = atom->nlocal + atom->nghost;
|
||||
const int nthreads = comm->nthreads;
|
||||
const int inum = list->inum;
|
||||
|
||||
// grow coordination array if necessary
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(NCo);
|
||||
nmax = atom->nmax;
|
||||
memory->create(NCo,nmax,"pair:NCo");
|
||||
}
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel default(shared)
|
||||
#endif
|
||||
{
|
||||
int ifrom, ito, tid;
|
||||
double **f;
|
||||
|
||||
f = loop_setup_thr(atom->f, ifrom, ito, tid, inum, nall, nthreads);
|
||||
|
||||
if (evflag) {
|
||||
if (eflag) {
|
||||
if (vflag_atom) eval<1,1,1>(f, ifrom, ito, tid);
|
||||
else eval<1,1,0>(f, ifrom, ito, tid);
|
||||
} else {
|
||||
if (vflag_atom) eval<1,0,1>(f, ifrom, ito, tid);
|
||||
else eval<1,0,0>(f, ifrom, ito, tid);
|
||||
}
|
||||
} else eval<0,0,0>(f, ifrom, ito, tid);
|
||||
|
||||
// reduce per thread forces into global force array.
|
||||
data_reduce_thr(&(atom->f[0][0]), nall, nthreads, 3, tid);
|
||||
} // end of omp parallel region
|
||||
|
||||
// reduce per thread energy and virial, if requested.
|
||||
if (evflag) ev_reduce_thr(this);
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
template <int EVFLAG, int EFLAG, int VFLAG_ATOM>
|
||||
void PairCombOMP::eval(double **f, int iifrom, int iito, int tid)
|
||||
{
|
||||
int i,j,k,ii,jj,kk,jnum,iparam_i;
|
||||
int itag,jtag,itype,jtype,ktype,iparam_ij,iparam_ijk;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double rsq,rsq1,rsq2;
|
||||
double delr1[3],delr2[3],fi[3],fj[3],fk[3];
|
||||
double zeta_ij,prefactor;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int mr1,mr2,mr3;
|
||||
int rsc,inty;
|
||||
double elp_ij,filp[3],fjlp[3],fklp[3];
|
||||
double iq,jq;
|
||||
double yaself;
|
||||
double potal,fac11,fac11e;
|
||||
double vionij,fvionij,sr1,sr2,sr3,Eov,Fov;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
|
||||
double **x = atom->x;
|
||||
double *q = atom->q;
|
||||
int *tag = atom->tag;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
yaself = vionij = fvionij = Eov = Fov = 0.0;
|
||||
|
||||
double fxtmp,fytmp,fztmp;
|
||||
double fjxtmp,fjytmp,fjztmp;
|
||||
|
||||
// self energy correction term: potal
|
||||
|
||||
potal_calc(potal,fac11,fac11e);
|
||||
|
||||
// loop over full neighbor list of my atoms
|
||||
|
||||
for (ii = iifrom; ii < iito; ++ii) {
|
||||
|
||||
i = ilist[ii];
|
||||
itag = tag[i];
|
||||
itype = map[type[i]];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
|
||||
iq = q[i];
|
||||
NCo[i] = 0;
|
||||
iparam_i = elem2param[itype][itype][itype];
|
||||
|
||||
// self energy, only on i atom
|
||||
|
||||
yaself = self(¶ms[iparam_i],iq,potal);
|
||||
|
||||
if (EVFLAG) ev_tally_thr(this,i,i,nlocal,0,yaself,
|
||||
0.0,0.0,0.0,0.0,0.0,tid);
|
||||
|
||||
// two-body interactions (long and short repulsive)
|
||||
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
jtag = tag[j];
|
||||
|
||||
if (itag > jtag) {
|
||||
if ((itag+jtag) % 2 == 0) continue;
|
||||
} else if (itag < jtag) {
|
||||
if ((itag+jtag) % 2 == 1) continue;
|
||||
} else {
|
||||
if (x[j][2] < ztmp) continue;
|
||||
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
|
||||
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
|
||||
}
|
||||
|
||||
// Qj calculates 2-body Coulombic
|
||||
|
||||
jtype = map[type[j]];
|
||||
jq = q[j];
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
|
||||
iparam_ij = elem2param[itype][jtype][jtype];
|
||||
|
||||
// long range q-dependent
|
||||
|
||||
if (rsq > params[iparam_ij].lcutsq) continue;
|
||||
|
||||
inty = intype[itype][jtype];
|
||||
|
||||
// polynomial three-point interpolation
|
||||
|
||||
tri_point(rsq, mr1, mr2, mr3, sr1, sr2, sr3, itype);
|
||||
|
||||
// 1/r energy and forces
|
||||
|
||||
direct(inty,mr1,mr2,mr3,rsq,sr1,sr2,sr3,iq,jq,
|
||||
potal,fac11,fac11e,vionij,fvionij);
|
||||
|
||||
// field correction to self energy
|
||||
|
||||
field(¶ms[iparam_ij],rsq,iq,jq,vionij,fvionij);
|
||||
|
||||
// polarization field
|
||||
// sums up long range forces
|
||||
|
||||
fxtmp += delx*fvionij;
|
||||
fytmp += dely*fvionij;
|
||||
fztmp += delz*fvionij;
|
||||
f[j][0] -= delx*fvionij;
|
||||
f[j][1] -= dely*fvionij;
|
||||
f[j][2] -= delz*fvionij;
|
||||
|
||||
if (EVFLAG)
|
||||
ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
|
||||
0.0,vionij,fvionij,delx,dely,delz,tid);
|
||||
|
||||
// short range q-independent
|
||||
|
||||
if (rsq > params[iparam_ij].cutsq) continue;
|
||||
|
||||
repulsive(¶ms[iparam_ij],rsq,fpair,EFLAG,evdwl,iq,jq);
|
||||
|
||||
// repulsion is pure two-body, sums up pair repulsive forces
|
||||
|
||||
fxtmp += delx*fpair;
|
||||
fytmp += dely*fpair;
|
||||
fztmp += delz*fpair;
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
|
||||
if (EVFLAG)
|
||||
ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
|
||||
evdwl,0.0,fpair,delx,dely,delz,tid);
|
||||
}
|
||||
|
||||
// accumulate coordination number information
|
||||
|
||||
if (cor_flag) {
|
||||
int numcoor = 0;
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
jtype = map[type[j]];
|
||||
iparam_ij = elem2param[itype][jtype][jtype];
|
||||
|
||||
if(params[iparam_ij].hfocor > 0.0 ) {
|
||||
delr1[0] = x[j][0] - xtmp;
|
||||
delr1[1] = x[j][1] - ytmp;
|
||||
delr1[2] = x[j][2] - ztmp;
|
||||
rsq1 = vec3_dot(delr1,delr1);
|
||||
|
||||
if (rsq1 > params[iparam_ij].cutsq) continue;
|
||||
++numcoor;
|
||||
}
|
||||
NCo[i] = numcoor;
|
||||
}
|
||||
}
|
||||
|
||||
// three-body interactions
|
||||
// skip immediately if I-J is not within cutoff
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
jtype = map[type[j]];
|
||||
iparam_ij = elem2param[itype][jtype][jtype];
|
||||
|
||||
// this Qj for q-dependent BSi
|
||||
|
||||
jq = q[j];
|
||||
|
||||
delr1[0] = x[j][0] - xtmp;
|
||||
delr1[1] = x[j][1] - ytmp;
|
||||
delr1[2] = x[j][2] - ztmp;
|
||||
rsq1 = vec3_dot(delr1,delr1);
|
||||
|
||||
if (rsq1 > params[iparam_ij].cutsq) continue;
|
||||
|
||||
// accumulate bondorder zeta for each i-j interaction via loop over k
|
||||
|
||||
fjxtmp = fjytmp = fjztmp = 0.0;
|
||||
zeta_ij = 0.0;
|
||||
cuo_flag1 = 0; cuo_flag2 = 0;
|
||||
|
||||
for (kk = 0; kk < jnum; kk++) {
|
||||
if (jj == kk) continue;
|
||||
k = jlist[kk];
|
||||
k &= NEIGHMASK;
|
||||
ktype = map[type[k]];
|
||||
iparam_ijk = elem2param[itype][jtype][ktype];
|
||||
|
||||
delr2[0] = x[k][0] - xtmp;
|
||||
delr2[1] = x[k][1] - ytmp;
|
||||
delr2[2] = x[k][2] - ztmp;
|
||||
rsq2 = vec3_dot(delr2,delr2);
|
||||
|
||||
if (rsq2 > params[iparam_ijk].cutsq) continue;
|
||||
|
||||
zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
|
||||
|
||||
if (params[iparam_ijk].hfocor == -2.0) cuo_flag1 = 1;
|
||||
if (params[iparam_ijk].hfocor == -1.0) cuo_flag2 = 1;
|
||||
}
|
||||
|
||||
if (cuo_flag1 && cuo_flag2) cuo_flag = 1;
|
||||
else cuo_flag = 0;
|
||||
|
||||
// pairwise force due to zeta
|
||||
|
||||
force_zeta(¶ms[iparam_ij],rsq1,zeta_ij,fpair,
|
||||
prefactor,EFLAG,evdwl,iq,jq);
|
||||
|
||||
// over-coordination correction for HfO2
|
||||
|
||||
if (cor_flag && NCo[i] != 0)
|
||||
Over_cor(¶ms[iparam_ij],rsq1,NCo[i],Eov, Fov);
|
||||
evdwl += Eov;
|
||||
fpair += Fov;
|
||||
|
||||
fxtmp += delr1[0]*fpair;
|
||||
fytmp += delr1[1]*fpair;
|
||||
fztmp += delr1[2]*fpair;
|
||||
fjxtmp -= delr1[0]*fpair;
|
||||
fjytmp -= delr1[1]*fpair;
|
||||
fjztmp -= delr1[2]*fpair;
|
||||
|
||||
if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,evdwl,0.0,
|
||||
-fpair,-delr1[0],-delr1[1],-delr1[2],tid);
|
||||
|
||||
// attractive term via loop over k (3-body forces)
|
||||
|
||||
for (kk = 0; kk < jnum; kk++) {
|
||||
if (jj == kk) continue;
|
||||
k = jlist[kk];
|
||||
k &= NEIGHMASK;
|
||||
ktype = map[type[k]];
|
||||
iparam_ijk = elem2param[itype][jtype][ktype];
|
||||
|
||||
delr2[0] = x[k][0] - xtmp;
|
||||
delr2[1] = x[k][1] - ytmp;
|
||||
delr2[2] = x[k][2] - ztmp;
|
||||
rsq2 = vec3_dot(delr2,delr2);
|
||||
if (rsq2 > params[iparam_ijk].cutsq) continue;
|
||||
|
||||
for (rsc = 0; rsc < 3; rsc++)
|
||||
fi[rsc] = fj[rsc] = fk[rsc] = 0.0;
|
||||
|
||||
attractive(¶ms[iparam_ijk],prefactor,
|
||||
rsq1,rsq2,delr1,delr2,fi,fj,fk);
|
||||
|
||||
// 3-body LP and BB correction and forces
|
||||
|
||||
elp_ij = elp(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
|
||||
flp(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2,filp,fjlp,fklp);
|
||||
|
||||
fxtmp += fi[0] + filp[0];
|
||||
fytmp += fi[1] + filp[1];
|
||||
fztmp += fi[2] + filp[2];
|
||||
fjxtmp += fj[0] + fjlp[0];
|
||||
fjytmp += fj[1] + fjlp[1];
|
||||
fjztmp += fj[2] + fjlp[2];
|
||||
f[k][0] += fk[0] + fklp[0];
|
||||
f[k][1] += fk[1] + fklp[1];
|
||||
f[k][2] += fk[2] + fklp[2];
|
||||
|
||||
if (EVFLAG)
|
||||
ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
|
||||
elp_ij,0.0,0.0,0.0,0.0,0.0, tid);
|
||||
if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,tid);
|
||||
}
|
||||
f[j][0] += fjxtmp;
|
||||
f[j][1] += fjytmp;
|
||||
f[j][2] += fjztmp;
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
|
||||
if (cuo_flag) params[iparam_i].cutsq *= 0.65;
|
||||
}
|
||||
cuo_flag = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairCombOMP::yasu_char(double *qf_fix, int &igroup)
|
||||
{
|
||||
int ii;
|
||||
double potal,fac11,fac11e;
|
||||
|
||||
const double * const * const x = atom->x;
|
||||
const double * const q = atom->q;
|
||||
const int * const type = atom->type;
|
||||
|
||||
const int inum = list->inum;
|
||||
const int * const ilist = list->ilist;
|
||||
const int * const numneigh = list->numneigh;
|
||||
const int * const * const firstneigh = list->firstneigh;
|
||||
|
||||
const int * const mask = atom->mask;
|
||||
const int groupbit = group->bitmask[igroup];
|
||||
|
||||
qf = qf_fix;
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
const int i = ilist[ii];
|
||||
if (mask[i] & groupbit)
|
||||
qf[i] = 0.0;
|
||||
}
|
||||
|
||||
// communicating charge force to all nodes, first forward then reverse
|
||||
|
||||
comm->forward_comm_pair(this);
|
||||
|
||||
// self energy correction term: potal
|
||||
|
||||
potal_calc(potal,fac11,fac11e);
|
||||
|
||||
// loop over full neighbor list of my atoms
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel for private(ii) default(none) shared(potal,fac11e)
|
||||
#endif
|
||||
for (ii = 0; ii < inum; ii ++) {
|
||||
double fqi,fqj,fqij,fqji,fqjj,delr1[3],delr2[3];
|
||||
double sr1,sr2,sr3;
|
||||
int mr1,mr2,mr3;
|
||||
|
||||
const int i = ilist[ii];
|
||||
|
||||
if (mask[i] & groupbit) {
|
||||
fqi = fqj = fqij = fqji = fqjj = 0.0; // should not be needed.
|
||||
int itype = map[type[i]];
|
||||
const double xtmp = x[i][0];
|
||||
const double ytmp = x[i][1];
|
||||
const double ztmp = x[i][2];
|
||||
const double iq = q[i];
|
||||
const int iparam_i = elem2param[itype][itype][itype];
|
||||
|
||||
// charge force from self energy
|
||||
|
||||
fqi = qfo_self(¶ms[iparam_i],iq,potal);
|
||||
|
||||
// two-body interactions
|
||||
|
||||
const int * const jlist = firstneigh[i];
|
||||
const int jnum = numneigh[i];
|
||||
|
||||
for (int jj = 0; jj < jnum; jj++) {
|
||||
const int j = jlist[jj] & NEIGHMASK;
|
||||
const int jtype = map[type[j]];
|
||||
double jq = q[j];
|
||||
|
||||
delr1[0] = x[j][0] - xtmp;
|
||||
delr1[1] = x[j][1] - ytmp;
|
||||
delr1[2] = x[j][2] - ztmp;
|
||||
double rsq1 = vec3_dot(delr1,delr1);
|
||||
|
||||
const int iparam_ij = elem2param[itype][jtype][jtype];
|
||||
|
||||
// long range q-dependent
|
||||
|
||||
if (rsq1 > params[iparam_ij].lcutsq) continue;
|
||||
|
||||
const int inty = intype[itype][jtype];
|
||||
|
||||
// polynomial three-point interpolation
|
||||
|
||||
tri_point(rsq1,mr1,mr2,mr3,sr1,sr2,sr3,itype);
|
||||
|
||||
// 1/r charge forces
|
||||
|
||||
qfo_direct(inty,mr1,mr2,mr3,rsq1,sr1,sr2,sr3,fac11e,fqij);
|
||||
|
||||
// field correction to self energy and charge force
|
||||
|
||||
qfo_field(¶ms[iparam_ij],rsq1,iq,jq,fqji,fqjj);
|
||||
fqi += jq * fqij + fqji;
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp atomic
|
||||
#endif
|
||||
qf[j] += (iq * fqij + fqjj);
|
||||
|
||||
// polarization field charge force
|
||||
// three-body interactions
|
||||
|
||||
if (rsq1 > params[iparam_ij].cutsq) continue;
|
||||
|
||||
double zeta_ij = 0.0;
|
||||
|
||||
for (int kk = 0; kk < jnum; kk++) {
|
||||
if (jj == kk) continue;
|
||||
const int k = jlist[kk] & NEIGHMASK;
|
||||
const int ktype = map[type[k]];
|
||||
const int iparam_ijk = elem2param[itype][jtype][ktype];
|
||||
|
||||
delr2[0] = x[k][0] - xtmp;
|
||||
delr2[1] = x[k][1] - ytmp;
|
||||
delr2[2] = x[k][2] - ztmp;
|
||||
const double rsq2 = vec3_dot(delr2,delr2);
|
||||
|
||||
if (rsq2 > params[iparam_ijk].cutsq) continue;
|
||||
zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
|
||||
}
|
||||
|
||||
// charge force in Aij and Bij
|
||||
|
||||
qfo_short(¶ms[iparam_ij],rsq1,zeta_ij,iq,jq,fqij,fqjj);
|
||||
fqi += fqij;
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp atomic
|
||||
#endif
|
||||
qf[j] += fqjj;
|
||||
}
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp atomic
|
||||
#endif
|
||||
qf[i] += fqi;
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
comm->reverse_comm_pair(this);
|
||||
|
||||
// sum charge force on each node and return it
|
||||
|
||||
double eneg = 0.0;
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
const int i = ilist[ii];
|
||||
if (mask[i] & groupbit)
|
||||
eneg += qf[i];
|
||||
}
|
||||
double enegtot;
|
||||
MPI_Allreduce(&eneg,&enegtot,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
return enegtot;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairCombOMP::memory_usage()
|
||||
{
|
||||
double bytes = memory_usage_thr();
|
||||
bytes += PairComb::memory_usage();
|
||||
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -1,45 +0,0 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(comb/omp,PairCombOMP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_COMB_OMP_H
|
||||
#define LMP_PAIR_COMB_OMP_H
|
||||
|
||||
#include "pair_comb.h"
|
||||
#include "thr_omp.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairCombOMP : public PairComb, public ThrOMP {
|
||||
|
||||
public:
|
||||
PairCombOMP(class LAMMPS *);
|
||||
|
||||
virtual void compute(int, int);
|
||||
virtual double memory_usage();
|
||||
|
||||
virtual double yasu_char(double *, int &);
|
||||
|
||||
private:
|
||||
template <int EVFLAG, int EFLAG, int VFLAG_ATOM>
|
||||
void eval(double **f, int ifrom, int ito, int tid);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
Loading…
Reference in New Issue