Page MenuHomec4science
Files F66033220

pair_snap.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Jun 7, 20:34

pair_snap.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "pair_snap.h"
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "sna.h"
#include "openmp_snap.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#include <cmath>
using namespace LAMMPS_NS;
#define MAXLINE 1024
#define MAXWORD 3
// Outstanding issues with quadratic term
// 1. there seems to a problem with compute_optimized energy calc
// it does not match compute_regular, even when quadratic coeffs = 0
/* ---------------------------------------------------------------------- */
PairSNAP::PairSNAP(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
nelements = 0;
elements = NULL;
radelem = NULL;
wjelem = NULL;
coeffelem = NULL;
nmax = 0;
nthreads = 1;
schedule_user = 0;
schedule_time_guided = -1;
schedule_time_dynamic = -1;
ncalls_neigh =-1;
ilistmask_max = 0;
ilistmask = NULL;
ghostinum = 0;
ghostilist_max = 0;
ghostilist = NULL;
ghostnumneigh_max = 0;
ghostnumneigh = NULL;
ghostneighs = NULL;
ghostfirstneigh = NULL;
ghostneighs_total = 0;
ghostneighs_max = 0;
i_max = 0;
i_neighmax = 0;
i_numpairs = 0;
i_rij = NULL;
i_inside = NULL;
i_wj = NULL;
i_rcutij = NULL;
i_ninside = NULL;
i_pairs = NULL;
i_uarraytot_r = NULL;
i_uarraytot_i = NULL;
i_zarray_r = NULL;
i_zarray_i =NULL;
use_shared_arrays = 0;
#ifdef TIMING_INFO
timers[0] = 0;
timers[1] = 0;
timers[2] = 0;
timers[3] = 0;
#endif
// Need to set this because restart not handled by PairHybrid
sna = NULL;
}
/* ---------------------------------------------------------------------- */
PairSNAP::~PairSNAP()
{
if (nelements) {
for (int i = 0; i < nelements; i++)
delete[] elements[i];
delete[] elements;
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(coeffelem);
}
// Need to set this because restart not handled by PairHybrid
if (sna) {
#ifdef TIMING_INFO
double time[5];
double timeave[5];
double timeave_mpi[5];
double timemax_mpi[5];
for (int i = 0; i < 5; i++) {
time[i] = 0;
timeave[i] = 0;
for (int tid = 0; tid<nthreads; tid++) {
if (sna[tid]->timers[i]>time[i])
time[i] = sna[tid]->timers[i];
timeave[i] += sna[tid]->timers[i];
}
timeave[i] /= nthreads;
}
MPI_Reduce(timeave, timeave_mpi, 5, MPI_DOUBLE, MPI_SUM, 0, world);
MPI_Reduce(time, timemax_mpi, 5, MPI_DOUBLE, MPI_MAX, 0, world);
#endif
for (int tid = 0; tid<nthreads; tid++)
delete sna[tid];
delete [] sna;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(map);
}
}
void PairSNAP::compute(int eflag, int vflag)
{
if (use_optimized)
compute_optimized(eflag, vflag);
else
compute_regular(eflag, vflag);
}
/* ----------------------------------------------------------------------
This version is a straightforward implementation
---------------------------------------------------------------------- */
void PairSNAP::compute_regular(int eflag, int vflag)
{
int i,j,jnum,ninside;
double delx,dely,delz,evdwl,rsq;
double fij[3];
int *jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
class SNA* snaptr = sna[0];
numneigh = list->numneigh;
firstneigh = list->firstneigh;
for (int ii = 0; ii < list->inum; ii++) {
i = list->ilist[ii];
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
const int itype = type[i];
const int ielem = map[itype];
const double radi = radelem[ielem];
jlist = firstneigh[i];
jnum = numneigh[i];
// insure rij, inside, wj, and rcutij are of size jnum
snaptr->grow_rij(jnum);
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
// wj = weights for neighbors of I within cutoff
// rcutij = cutoffs for neighbors of I within cutoff
// note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = x[j][0] - xtmp;
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j];
int jelem = map[jtype];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
snaptr->rij[ninside][0] = delx;
snaptr->rij[ninside][1] = dely;
snaptr->rij[ninside][2] = delz;
snaptr->inside[ninside] = j;
snaptr->wj[ninside] = wjelem[jelem];
snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac;
ninside++;
}
}
// compute Ui, Zi, and Bi for atom I
snaptr->compute_ui(ninside);
snaptr->compute_zi();
// for neighbors of I within cutoff:
// compute dUi/drj and dBi/drj
// Fij = dEi/dRj = -dEi/dRi => add to Fi, subtract from Fj
double* coeffi = coeffelem[ielem];
for (int jj = 0; jj < ninside; jj++) {
int j = snaptr->inside[jj];
snaptr->compute_duidrj(snaptr->rij[jj],
snaptr->wj[jj],snaptr->rcutij[jj]);
snaptr->compute_dbidrj();
snaptr->copy_dbi2dbvec();
fij[0] = 0.0;
fij[1] = 0.0;
fij[2] = 0.0;
// linear contributions
for (int k = 1; k <= ncoeff; k++) {
double bgb = coeffi[k];
fij[0] += bgb*snaptr->dbvec[k-1][0];
fij[1] += bgb*snaptr->dbvec[k-1][1];
fij[2] += bgb*snaptr->dbvec[k-1][2];
}
// quadratic contributions
if (quadraticflag) {
snaptr->compute_bi();
snaptr->copy_bi2bvec();
int k = ncoeff+1;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bveci = snaptr->bvec[icoeff];
double fack = coeffi[k]*bveci;
double* dbveci = snaptr->dbvec[icoeff];
fij[0] += fack*snaptr->dbvec[icoeff][0];
fij[1] += fack*snaptr->dbvec[icoeff][1];
fij[2] += fack*snaptr->dbvec[icoeff][2];
k++;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
double facki = coeffi[k]*bveci;
double fackj = coeffi[k]*snaptr->bvec[jcoeff];
double* dbvecj = snaptr->dbvec[jcoeff];
fij[0] += facki*dbvecj[0]+fackj*dbveci[0];
fij[1] += facki*dbvecj[1]+fackj*dbveci[1];
fij[2] += facki*dbvecj[2]+fackj*dbveci[2];
k++;
}
}
}
f[i][0] += fij[0];
f[i][1] += fij[1];
f[i][2] += fij[2];
f[j][0] -= fij[0];
f[j][1] -= fij[1];
f[j][2] -= fij[2];
// tally per-atom virial contribution
if (vflag)
ev_tally_xyz(i,j,nlocal,newton_pair,0.0,0.0,
fij[0],fij[1],fij[2],
-snaptr->rij[jj][0],-snaptr->rij[jj][1],
-snaptr->rij[jj][2]);
}
// tally energy contribution
if (eflag) {
// evdwl = energy of atom I, sum over coeffs_k * Bi_k
evdwl = coeffi[0];
if (!quadraticflag) {
snaptr->compute_bi();
snaptr->copy_bi2bvec();
}
// E = beta.B + 0.5*B^t.alpha.B
// coeff[k] = beta[k-1] or
// coeff[k] = alpha_ii or
// coeff[k] = alpha_ij = alpha_ji, j != i
// linear contributions
for (int k = 1; k <= ncoeff; k++)
evdwl += coeffi[k]*snaptr->bvec[k-1];
// quadratic contributions
if (quadraticflag) {
int k = ncoeff+1;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bveci = snaptr->bvec[icoeff];
evdwl += 0.5*coeffi[k++]*bveci*bveci;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
evdwl += coeffi[k++]*bveci*snaptr->bvec[jcoeff];
}
}
}
ev_tally_full(i,2.0*evdwl,0.0,0.0,0.0,0.0,0.0);
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
This version is optimized for threading, micro-load balancing
---------------------------------------------------------------------- */
void PairSNAP::compute_optimized(int eflag, int vflag)
{
// if reneighboring took place do load_balance if requested
if (do_load_balance > 0 &&
(neighbor->ncalls != ncalls_neigh)) {
ghostinum = 0;
// reset local ghost neighbor lists
ncalls_neigh = neighbor->ncalls;
if (ilistmask_max < list->inum) {
memory->grow(ilistmask,list->inum,"PairSnap::ilistmask");
ilistmask_max = list->inum;
}
for (int i = 0; i < list->inum; i++)
ilistmask[i] = 1;
//multiple passes for loadbalancing
for (int i = 0; i < do_load_balance; i++)
load_balance();
}
int numpairs = 0;
for (int ii = 0; ii < list->inum; ii++) {
if ((do_load_balance <= 0) || ilistmask[ii]) {
int i = list->ilist[ii];
int jnum = list->numneigh[i];
numpairs += jnum;
}
}
if (do_load_balance)
for (int ii = 0; ii < ghostinum; ii++) {
int i = ghostilist[ii];
int jnum = ghostnumneigh[i];
numpairs += jnum;
}
// optimized schedule setting
int time_dynamic = 0;
int time_guided = 0;
if (schedule_user == 0) schedule_user = 4;
switch (schedule_user) {
case 1:
omp_set_schedule(omp_sched_static,1);
break;
case 2:
omp_set_schedule(omp_sched_dynamic,1);
break;
case 3:
omp_set_schedule(omp_sched_guided,2);
break;
case 4:
omp_set_schedule(omp_sched_auto,0);
break;
case 5:
if (numpairs < 8*nthreads) omp_set_schedule(omp_sched_dynamic,1);
else if (schedule_time_guided < 0.0) {
omp_set_schedule(omp_sched_guided,2);
if (!eflag && !vflag) time_guided = 1;
} else if (schedule_time_dynamic<0.0) {
omp_set_schedule(omp_sched_dynamic,1);
if (!eflag && !vflag) time_dynamic = 1;
} else if (schedule_time_guided<schedule_time_dynamic)
omp_set_schedule(omp_sched_guided,2);
else
omp_set_schedule(omp_sched_dynamic,1);
break;
}
if (use_shared_arrays)
build_per_atom_arrays();
#if defined(_OPENMP)
#pragma omp parallel shared(eflag,vflag,time_dynamic,time_guided) firstprivate(numpairs) default(none)
#endif
{
// begin of pragma omp parallel
int tid = omp_get_thread_num();
int** pairs_tid_unique = NULL;
int** pairs;
if (use_shared_arrays) pairs = i_pairs;
else {
memory->create(pairs_tid_unique,numpairs,4,"numpairs");
pairs = pairs_tid_unique;
}
if (!use_shared_arrays) {
numpairs = 0;
for (int ii = 0; ii < list->inum; ii++) {
if ((do_load_balance <= 0) || ilistmask[ii]) {
int i = list->ilist[ii];
int jnum = list->numneigh[i];
for (int jj = 0; jj<jnum; jj++) {
pairs[numpairs][0] = i;
pairs[numpairs][1] = jj;
pairs[numpairs][2] = -1;
numpairs++;
}
}
}
for (int ii = 0; ii < ghostinum; ii++) {
int i = ghostilist[ii];
int jnum = ghostnumneigh[i];
for (int jj = 0; jj<jnum; jj++) {
pairs[numpairs][0] = i;
pairs[numpairs][1] = jj;
pairs[numpairs][2] = -1;
numpairs++;
}
}
}
int ielem;
int jj,k,jnum,jtype,ninside;
double delx,dely,delz,evdwl,rsq;
double fij[3];
int *jlist,*numneigh,**firstneigh;
evdwl = 0.0;
#if defined(_OPENMP)
#pragma omp master
#endif
{
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
}
#if defined(_OPENMP)
#pragma omp barrier
{ ; }
#endif
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
#ifdef TIMING_INFO
// only update micro timers after setup
static int count=0;
if (count<2) {
sna[tid]->timers[0] = 0;
sna[tid]->timers[1] = 0;
sna[tid]->timers[2] = 0;
sna[tid]->timers[3] = 0;
sna[tid]->timers[4] = 0;
}
count++;
#endif
// did thread start working on interactions of new atom
int iold = -1;
double starttime, endtime;
if (time_dynamic || time_guided)
starttime = MPI_Wtime();
#if defined(_OPENMP)
#pragma omp for schedule(runtime)
#endif
for (int iijj = 0; iijj < numpairs; iijj++) {
int i = 0;
if (use_shared_arrays) {
i = i_pairs[iijj][0];
if (iold != i) {
set_sna_to_shared(tid,i_pairs[iijj][3]);
ielem = map[type[i]];
}
iold = i;
} else {
i = pairs[iijj][0];
if (iold != i) {
iold = i;
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
const int itype = type[i];
ielem = map[itype];
const double radi = radelem[ielem];
if (i < nlocal) {
jlist = firstneigh[i];
jnum = numneigh[i];
} else {
jlist = ghostneighs+ghostfirstneigh[i];
jnum = ghostnumneigh[i];
}
// insure rij, inside, wj, and rcutij are of size jnum
sna[tid]->grow_rij(jnum);
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
// wj = weights of neighbors of I within cutoff
// rcutij = cutoffs of neighbors of I within cutoff
// note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
ninside = 0;
for (jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
delx = x[j][0] - xtmp; //unitialised
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
int jelem = map[jtype];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) { //unitialised
sna[tid]->rij[ninside][0] = delx;
sna[tid]->rij[ninside][1] = dely;
sna[tid]->rij[ninside][2] = delz;
sna[tid]->inside[ninside] = j;
sna[tid]->wj[ninside] = wjelem[jelem];
sna[tid]->rcutij[ninside] = (radi + radelem[jelem])*rcutfac;
ninside++;
// update index list with inside index
pairs[iijj + (jj - pairs[iijj][1])][2] =
ninside-1; //unitialised
}
}
// compute Ui and Zi for atom I
sna[tid]->compute_ui(ninside); //unitialised
sna[tid]->compute_zi();
}
}
// for neighbors of I within cutoff:
// compute dUi/drj and dBi/drj
// Fij = dEi/dRj = -dEi/dRi => add to Fi, subtract from Fj
// entry into loop if inside index is set
double* coeffi = coeffelem[ielem];
if (pairs[iijj][2] >= 0) {
jj = pairs[iijj][2];
int j = sna[tid]->inside[jj];
sna[tid]->compute_duidrj(sna[tid]->rij[jj],
sna[tid]->wj[jj],sna[tid]->rcutij[jj]);
sna[tid]->compute_dbidrj();
sna[tid]->copy_dbi2dbvec();
fij[0] = 0.0;
fij[1] = 0.0;
fij[2] = 0.0;
// linear contributions
for (k = 1; k <= ncoeff; k++) {
double bgb = coeffi[k];
fij[0] += bgb*sna[tid]->dbvec[k-1][0];
fij[1] += bgb*sna[tid]->dbvec[k-1][1];
fij[2] += bgb*sna[tid]->dbvec[k-1][2];
}
// quadratic contributions
if (quadraticflag) {
sna[tid]->compute_bi();
sna[tid]->copy_bi2bvec();
int k = ncoeff+1;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bveci = sna[tid]->bvec[icoeff];
double fack = coeffi[k]*bveci;
double* dbveci = sna[tid]->dbvec[icoeff];
fij[0] += fack*sna[tid]->dbvec[icoeff][0];
fij[1] += fack*sna[tid]->dbvec[icoeff][1];
fij[2] += fack*sna[tid]->dbvec[icoeff][2];
k++;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
double facki = coeffi[k]*bveci;
double fackj = coeffi[k]*sna[tid]->bvec[jcoeff];
double* dbvecj = sna[tid]->dbvec[jcoeff];
fij[0] += facki*dbvecj[0]+fackj*dbveci[0];
fij[1] += facki*dbvecj[1]+fackj*dbveci[1];
fij[2] += facki*dbvecj[2]+fackj*dbveci[2];
k++;
}
}
}
#if defined(_OPENMP)
#pragma omp critical
#endif
{
f[i][0] += fij[0];
f[i][1] += fij[1];
f[i][2] += fij[2];
f[j][0] -= fij[0];
f[j][1] -= fij[1];
f[j][2] -= fij[2];
// tally per-atom virial contribution
if (vflag)
ev_tally_xyz(i,j,nlocal,newton_pair,0.0,0.0,
fij[0],fij[1],fij[2],
-sna[tid]->rij[jj][0],-sna[tid]->rij[jj][1],
-sna[tid]->rij[jj][2]);
}
}
// evdwl = energy of atom I, sum over coeffs_k * Bi_k
// only call this for first pair of each atom i
// if atom has no pairs, eatom=0, which is wrong
if (eflag&&pairs[iijj][1] == 0) {
evdwl = coeffi[0];
sna[tid]->compute_bi();
sna[tid]->copy_bi2bvec();
// E = beta.B + 0.5*B^t.alpha.B
// coeff[k] = beta[k-1] or
// coeff[k] = alpha_ii or
// coeff[k] = alpha_ij = alpha_ji, j != i
// linear contributions
for (int k = 1; k <= ncoeff; k++)
evdwl += coeffi[k]*sna[tid]->bvec[k-1];
// quadratic contributions
if (quadraticflag) {
int k = ncoeff+1;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bveci = sna[tid]->bvec[icoeff];
evdwl += 0.5*coeffi[k++]*bveci*bveci;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
evdwl += coeffi[k++]*bveci*sna[tid]->bvec[jcoeff];
}
}
}
#if defined(_OPENMP)
#pragma omp critical
#endif
ev_tally_full(i,2.0*evdwl,0.0,0.0,0.0,0.0,0.0);
}
}
if (time_dynamic || time_guided)
endtime = MPI_Wtime();
if (time_dynamic) schedule_time_dynamic = endtime - starttime;
if (time_guided) schedule_time_guided = endtime - starttime;
if (!use_shared_arrays) memory->destroy(pairs);
}// end of pragma omp parallel
if (vflag_fdotr) virial_fdotr_compute();
}
inline int PairSNAP::equal(double* x,double* y)
{
double dist2 =
(x[0]-y[0])*(x[0]-y[0]) +
(x[1]-y[1])*(x[1]-y[1]) +
(x[2]-y[2])*(x[2]-y[2]);
if (dist2 < 1e-20) return 1;
return 0;
}
inline double PairSNAP::dist2(double* x,double* y)
{
return
(x[0]-y[0])*(x[0]-y[0]) +
(x[1]-y[1])*(x[1]-y[1]) +
(x[2]-y[2])*(x[2]-y[2]);
}
// return extra communication cutoff
// extra_cutoff = max(subdomain_length)
double PairSNAP::extra_cutoff()
{
double sublo[3],subhi[3];
if (domain->triclinic == 0) {
for (int dim = 0 ; dim < 3 ; dim++) {
sublo[dim] = domain->sublo[dim];
subhi[dim] = domain->subhi[dim];
}
} else {
domain->lamda2x(domain->sublo_lamda,sublo);
domain->lamda2x(domain->subhi_lamda,subhi);
}
double sub_size[3];
for (int dim = 0; dim < 3; dim++)
sub_size[dim] = subhi[dim] - sublo[dim];
double max_sub_size = 0;
for (int dim = 0; dim < 3; dim++)
max_sub_size = MAX(max_sub_size,sub_size[dim]);
// note: for triclinic, probably need something different
// see Comm::setup()
return max_sub_size;
}
// micro load_balancer: each MPI process will
// check with each of its 26 neighbors,
// whether an imbalance exists in the number
// of atoms to calculate forces for.
// If it does it will set ilistmask of one of
// its local atoms to zero, and send its Tag
// to the neighbor process. The neighboring process
// will check its ghost list for the
// ghost atom with the same Tag which is closest
// to its domain center, and build a
// neighborlist for this ghost atom. For this to work,
// the communication cutoff has to be
// as large as the neighbor cutoff +
// maximum subdomain length.
// Note that at most one atom is exchanged per processor pair.
// Also note that the local atom assignment
// doesn't change. This load balancer will cause
// some ghost atoms to have full neighborlists
// which are unique to PairSNAP.
// They are not part of the generally accessible neighborlist.
// At the same time corresponding local atoms on
// other MPI processes will not be
// included in the force computation since
// their ilistmask is 0. This does not effect
// any other classes which might
// access the same general neighborlist.
// Reverse communication (newton on) of forces is required.
// Currently the load balancer does two passes,
// since its exchanging atoms upstream and downstream.
void PairSNAP::load_balance()
{
double sublo[3],subhi[3];
if (domain->triclinic == 0) {
double* sublotmp = domain->sublo;
double* subhitmp = domain->subhi;
for (int dim = 0 ; dim<3 ; dim++) {
sublo[dim]=sublotmp[dim];
subhi[dim]=subhitmp[dim];
}
} else {
double* sublotmp = domain->sublo_lamda;
double* subhitmp = domain->subhi_lamda;
domain->lamda2x(sublotmp,sublo);
domain->lamda2x(subhitmp,subhi);
}
//if (list->inum==0) list->grow(atom->nmax);
int nlocal = ghostinum;
for (int i=0; i < list->inum; i++)
if (ilistmask[i]) nlocal++;
int ***grid2proc = comm->grid2proc;
int* procgrid = comm->procgrid;
int nlocal_up,nlocal_down;
MPI_Request request;
double sub_mid[3];
for (int dim=0; dim<3; dim++)
sub_mid[dim] = (subhi[dim] + sublo[dim])/2;
if (comm->cutghostuser <
neighbor->cutneighmax+extra_cutoff())
error->all(FLERR,"Communication cutoff too small for SNAP micro load balancing");
int nrecv = ghostinum;
int totalsend = 0;
int nsend = 0;
int depth = 1;
for (int dx = -depth; dx < depth+1; dx++)
for (int dy = -depth; dy < depth+1; dy++)
for (int dz = -depth; dz < depth+1; dz++) {
if (dx == dy && dy == dz && dz == 0) continue;
int sendloc[3] = {comm->myloc[0],
comm->myloc[1], comm->myloc[2]
};
sendloc[0] += dx;
sendloc[1] += dy;
sendloc[2] += dz;
for (int dim = 0; dim < 3; dim++)
if (sendloc[dim] >= procgrid[dim])
sendloc[dim] = sendloc[dim] - procgrid[dim];
for (int dim = 0; dim < 3; dim++)
if (sendloc[dim] < 0)
sendloc[dim] = procgrid[dim] + sendloc[dim];
int recvloc[3] = {comm->myloc[0],
comm->myloc[1], comm->myloc[2]
};
recvloc[0] -= dx;
recvloc[1] -= dy;
recvloc[2] -= dz;
for (int dim = 0; dim < 3; dim++)
if (recvloc[dim] < 0)
recvloc[dim] = procgrid[dim] + recvloc[dim];
for (int dim = 0; dim < 3; dim++)
if (recvloc[dim] >= procgrid[dim])
recvloc[dim] = recvloc[dim] - procgrid[dim];
int sendproc = grid2proc[sendloc[0]][sendloc[1]][sendloc[2]];
int recvproc = grid2proc[recvloc[0]][recvloc[1]][recvloc[2]];
// two stage process, first upstream movement, then downstream
MPI_Sendrecv(&nlocal,1,MPI_INT,sendproc,0,
&nlocal_up,1,MPI_INT,recvproc,0,world,MPI_STATUS_IGNORE);
MPI_Sendrecv(&nlocal,1,MPI_INT,recvproc,0,
&nlocal_down,1,MPI_INT,sendproc,0,world,MPI_STATUS_IGNORE);
nsend = 0;
// send upstream
if (nlocal > nlocal_up+1) {
int i = totalsend++;
while(i < list->inum && ilistmask[i] == 0)
i = totalsend++;
if (i < list->inum)
MPI_Isend(&atom->tag[i],1,MPI_INT,recvproc,0,world,&request);
else {
int j = -1;
MPI_Isend(&j,1,MPI_INT,recvproc,0,world,&request);
}
if (i < list->inum) {
for (int j = 0; j < list->inum; j++)
if (list->ilist[j] == i)
ilistmask[j] = 0;
nsend = 1;
}
}
// recv downstream
if (nlocal < nlocal_down-1) {
nlocal++;
int get_tag = -1;
MPI_Recv(&get_tag,1,MPI_INT,sendproc,0,world,MPI_STATUS_IGNORE);
// if get_tag -1 the other process didnt have local atoms to send
if (get_tag >= 0) {
if (ghostinum >= ghostilist_max) {
memory->grow(ghostilist,ghostinum+10,
"PairSnap::ghostilist");
ghostilist_max = ghostinum+10;
}
if (atom->nlocal + atom->nghost >= ghostnumneigh_max) {
ghostnumneigh_max = atom->nlocal+atom->nghost+100;
memory->grow(ghostnumneigh,ghostnumneigh_max,
"PairSnap::ghostnumneigh");
memory->grow(ghostfirstneigh,ghostnumneigh_max,
"PairSnap::ghostfirstneigh");
}
// find closest ghost image of the transfered particle
double mindist = 1e200;
int closestghost = -1;
for (int j = 0; j < atom->nlocal + atom->nghost; j++)
if (atom->tag[j] == get_tag)
if (dist2(sub_mid, atom->x[j]) < mindist) {
closestghost = j;
mindist = dist2(sub_mid, atom->x[j]);
}
// build neighborlist for this particular
// ghost atom, and add it to list->ilist
if (ghostneighs_max - ghostneighs_total <
neighbor->oneatom) {
memory->grow(ghostneighs,
ghostneighs_total + neighbor->oneatom,
"PairSnap::ghostneighs");
ghostneighs_max = ghostneighs_total + neighbor->oneatom;
}
int j = closestghost;
ghostilist[ghostinum] = j;
ghostnumneigh[j] = 0;
ghostfirstneigh[j] = ghostneighs_total;
ghostinum++;
int* jlist = ghostneighs + ghostfirstneigh[j];
// find all neighbors by looping
// over all local and ghost atoms
for (int k = 0; k < atom->nlocal + atom->nghost; k++)
if (dist2(atom->x[j],atom->x[k]) <
neighbor->cutneighmax*neighbor->cutneighmax) {
jlist[ghostnumneigh[j]] = k;
ghostnumneigh[j]++;
ghostneighs_total++;
}
}
if (get_tag >= 0) nrecv++;
}
// decrease nlocal later, so that it is the
// initial number both for receiving and sending
if (nsend) nlocal--;
// second pass through the grid
MPI_Sendrecv(&nlocal,1,MPI_INT,sendproc,0,
&nlocal_up,1,MPI_INT,recvproc,0,world,MPI_STATUS_IGNORE);
MPI_Sendrecv(&nlocal,1,MPI_INT,recvproc,0,
&nlocal_down,1,MPI_INT,sendproc,0,world,MPI_STATUS_IGNORE);
// send downstream
nsend=0;
if (nlocal > nlocal_down+1) {
int i = totalsend++;
while(i < list->inum && ilistmask[i]==0) i = totalsend++;
if (i < list->inum)
MPI_Isend(&atom->tag[i],1,MPI_INT,sendproc,0,world,&request);
else {
int j =- 1;
MPI_Isend(&j,1,MPI_INT,sendproc,0,world,&request);
}
if (i < list->inum) {
for (int j=0; j<list->inum; j++)
if (list->ilist[j] == i) ilistmask[j] = 0;
nsend = 1;
}
}
// receive upstream
if (nlocal < nlocal_up-1) {
nlocal++;
int get_tag = -1;
MPI_Recv(&get_tag,1,MPI_INT,recvproc,0,world,MPI_STATUS_IGNORE);
if (get_tag >= 0) {
if (ghostinum >= ghostilist_max) {
memory->grow(ghostilist,ghostinum+10,
"PairSnap::ghostilist");
ghostilist_max = ghostinum+10;
}
if (atom->nlocal + atom->nghost >= ghostnumneigh_max) {
ghostnumneigh_max = atom->nlocal + atom->nghost + 100;
memory->grow(ghostnumneigh,ghostnumneigh_max,
"PairSnap::ghostnumneigh");
memory->grow(ghostfirstneigh,ghostnumneigh_max,
"PairSnap::ghostfirstneigh");
}
// find closest ghost image of the transfered particle
double mindist = 1e200;
int closestghost = -1;
for (int j = 0; j < atom->nlocal + atom->nghost; j++)
if (atom->tag[j] == get_tag)
if (dist2(sub_mid,atom->x[j])<mindist) {
closestghost = j;
mindist = dist2(sub_mid,atom->x[j]);
}
// build neighborlist for this particular ghost atom
if (ghostneighs_max-ghostneighs_total < neighbor->oneatom) {
memory->grow(ghostneighs,ghostneighs_total + neighbor->oneatom,
"PairSnap::ghostneighs");
ghostneighs_max = ghostneighs_total + neighbor->oneatom;
}
int j = closestghost;
ghostilist[ghostinum] = j;
ghostnumneigh[j] = 0;
ghostfirstneigh[j] = ghostneighs_total;
ghostinum++;
int* jlist = ghostneighs + ghostfirstneigh[j];
for (int k = 0; k < atom->nlocal + atom->nghost; k++)
if (dist2(atom->x[j],atom->x[k]) <
neighbor->cutneighmax*neighbor->cutneighmax) {
jlist[ghostnumneigh[j]] = k;
ghostnumneigh[j]++;
ghostneighs_total++;
}
}
if (get_tag >= 0) nrecv++;
}
if (nsend) nlocal--;
}
}
void PairSNAP::set_sna_to_shared(int snaid,int i)
{
sna[snaid]->rij = i_rij[i];
sna[snaid]->inside = i_inside[i];
sna[snaid]->wj = i_wj[i];
sna[snaid]->rcutij = i_rcutij[i];
sna[snaid]->zarray_r = i_zarray_r[i];
sna[snaid]->zarray_i = i_zarray_i[i];
sna[snaid]->uarraytot_r = i_uarraytot_r[i];
sna[snaid]->uarraytot_i = i_uarraytot_i[i];
}
void PairSNAP::build_per_atom_arrays()
{
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&starttime);
#endif
int count = 0;
int neighmax = 0;
for (int ii = 0; ii < list->inum; ii++)
if ((do_load_balance <= 0) || ilistmask[ii]) {
neighmax=MAX(neighmax,list->numneigh[list->ilist[ii]]);
++count;
}
for (int ii = 0; ii < ghostinum; ii++) {
neighmax=MAX(neighmax,ghostnumneigh[ghostilist[ii]]);
++count;
}
if (i_max < count || i_neighmax < neighmax) {
int i_maxt = MAX(count,i_max);
i_neighmax = MAX(neighmax,i_neighmax);
memory->destroy(i_rij);
memory->destroy(i_inside);
memory->destroy(i_wj);
memory->destroy(i_rcutij);
memory->destroy(i_ninside);
memory->destroy(i_pairs);
memory->create(i_rij,i_maxt,i_neighmax,3,"PairSNAP::i_rij");
memory->create(i_inside,i_maxt,i_neighmax,"PairSNAP::i_inside");
memory->create(i_wj,i_maxt,i_neighmax,"PairSNAP::i_wj");
memory->create(i_rcutij,i_maxt,i_neighmax,"PairSNAP::i_rcutij");
memory->create(i_ninside,i_maxt,"PairSNAP::i_ninside");
memory->create(i_pairs,i_maxt*i_neighmax,4,"PairSNAP::i_pairs");
}
if (i_max < count) {
int jdim = sna[0]->twojmax+1;
memory->destroy(i_uarraytot_r);
memory->destroy(i_uarraytot_i);
memory->create(i_uarraytot_r,count,jdim,jdim,jdim,
"PairSNAP::i_uarraytot_r");
memory->create(i_uarraytot_i,count,jdim,jdim,jdim,
"PairSNAP::i_uarraytot_i");
if (i_zarray_r != NULL)
for (int i = 0; i < i_max; i++) {
memory->destroy(i_zarray_r[i]);
memory->destroy(i_zarray_i[i]);
}
delete [] i_zarray_r;
delete [] i_zarray_i;
i_zarray_r = new double*****[count];
i_zarray_i = new double*****[count];
for (int i = 0; i < count; i++) {
memory->create(i_zarray_r[i],jdim,jdim,jdim,jdim,jdim,
"PairSNAP::i_zarray_r");
memory->create(i_zarray_i[i],jdim,jdim,jdim,jdim,jdim,
"PairSNAP::i_zarray_i");
}
}
if (i_max < count)
i_max = count;
count = 0;
i_numpairs = 0;
for (int ii = 0; ii < list->inum; ii++) {
if ((do_load_balance <= 0) || ilistmask[ii]) {
int i = list->ilist[ii];
int jnum = list->numneigh[i];
int* jlist = list->firstneigh[i];
const double xtmp = atom->x[i][0];
const double ytmp = atom->x[i][1];
const double ztmp = atom->x[i][2];
const int itype = atom->type[i];
const int ielem = map[itype];
const double radi = radelem[ielem];
int ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double delx = atom->x[j][0] - xtmp;
const double dely = atom->x[j][1] - ytmp;
const double delz = atom->x[j][2] - ztmp;
const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = atom->type[j];
int jelem = map[jtype];
i_pairs[i_numpairs][0] = i;
i_pairs[i_numpairs][1] = jj;
i_pairs[i_numpairs][2] = -1;
i_pairs[i_numpairs][3] = count;
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
i_rij[count][ninside][0] = delx;
i_rij[count][ninside][1] = dely;
i_rij[count][ninside][2] = delz;
i_inside[count][ninside] = j;
i_wj[count][ninside] = wjelem[jelem];
i_rcutij[count][ninside] = (radi + radelem[jelem])*rcutfac;
// update index list with inside index
i_pairs[i_numpairs][2] = ninside++;
}
i_numpairs++;
}
i_ninside[count] = ninside;
count++;
}
}
for (int ii = 0; ii < ghostinum; ii++) {
int i = ghostilist[ii];
int jnum = ghostnumneigh[i];
int* jlist = ghostneighs+ghostfirstneigh[i];
const double xtmp = atom->x[i][0];
const double ytmp = atom->x[i][1];
const double ztmp = atom->x[i][2];
const int itype = atom->type[i];
const int ielem = map[itype];
const double radi = radelem[ielem];
int ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double delx = atom->x[j][0] - xtmp;
const double dely = atom->x[j][1] - ytmp;
const double delz = atom->x[j][2] - ztmp;
const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = atom->type[j];
int jelem = map[jtype];
i_pairs[i_numpairs][0] = i;
i_pairs[i_numpairs][1] = jj;
i_pairs[i_numpairs][2] = -1;
i_pairs[i_numpairs][3] = count;
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
i_rij[count][ninside][0] = delx;
i_rij[count][ninside][1] = dely;
i_rij[count][ninside][2] = delz;
i_inside[count][ninside] = j;
i_wj[count][ninside] = wjelem[jelem];
i_rcutij[count][ninside] = (radi + radelem[jelem])*rcutfac;
// update index list with inside index
i_pairs[i_numpairs][2] = ninside++;
}
i_numpairs++;
}
i_ninside[count] = ninside;
count++;
}
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&endtime);
timers[0]+=(endtime.tv_sec-starttime.tv_sec+1.0*
(endtime.tv_nsec-starttime.tv_nsec)/1000000000);
#endif
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&starttime);
#endif
#if defined(_OPENMP)
#pragma omp parallel for shared(count) default(none)
#endif
for (int ii=0; ii < count; ii++) {
int tid = omp_get_thread_num();
set_sna_to_shared(tid,ii);
//sna[tid]->compute_ui(i_ninside[ii]);
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&starttime);
#endif
sna[tid]->compute_ui_omp(i_ninside[ii],MAX(int(nthreads/count),1));
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&endtime);
sna[tid]->timers[0]+=(endtime.tv_sec-starttime.tv_sec+1.0*
(endtime.tv_nsec-starttime.tv_nsec)/1000000000);
#endif
}
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&starttime);
#endif
for (int ii=0; ii < count; ii++) {
int tid = 0;//omp_get_thread_num();
set_sna_to_shared(tid,ii);
sna[tid]->compute_zi_omp(MAX(int(nthreads/count),1));
}
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&endtime);
sna[0]->timers[1]+=(endtime.tv_sec-starttime.tv_sec+1.0*
(endtime.tv_nsec-starttime.tv_nsec)/1000000000);
#endif
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME,&endtime);
timers[1]+=(endtime.tv_sec-starttime.tv_sec+1.0*
(endtime.tv_nsec-starttime.tv_nsec)/1000000000);
#endif
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSNAP::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(map,n+1,"pair:map");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSNAP::settings(int narg, char **arg)
{
// set default values for optional arguments
nthreads = -1;
use_shared_arrays=-1;
do_load_balance = 0;
use_optimized = 1;
// optional arguments
for (int i=0; i < narg; i++) {
if (i+2>narg) error->all(FLERR,"Illegal pair_style command");
if (strcmp(arg[i],"nthreads")==0) {
nthreads=force->inumeric(FLERR,arg[++i]);
#if defined(LMP_USER_OMP)
error->all(FLERR,"Must set number of threads via package omp command");
#else
omp_set_num_threads(nthreads);
comm->nthreads=nthreads;
#endif
continue;
}
if (strcmp(arg[i],"optimized")==0) {
use_optimized=force->inumeric(FLERR,arg[++i]);
continue;
}
if (strcmp(arg[i],"shared")==0) {
use_shared_arrays=force->inumeric(FLERR,arg[++i]);
continue;
}
if (strcmp(arg[i],"loadbalance")==0) {
do_load_balance = force->inumeric(FLERR,arg[++i]);
if (do_load_balance) {
double mincutoff = extra_cutoff() +
rcutmax + neighbor->skin;
if (comm->cutghostuser < mincutoff) {
char buffer[255];
//apparently mincutoff is 0 after sprintf command ?????
double tmp = mincutoff + 0.1;
sprintf(buffer, "Communication cutoff is too small "
"for SNAP micro load balancing, increased to %lf",
mincutoff+0.1);
if (comm->me==0)
error->warning(FLERR,buffer);
comm->cutghostuser = tmp;
}
}
continue;
}
if (strcmp(arg[i],"schedule")==0) {
i++;
if (strcmp(arg[i],"static")==0)
schedule_user = 1;
if (strcmp(arg[i],"dynamic")==0)
schedule_user = 2;
if (strcmp(arg[i],"guided")==0)
schedule_user = 3;
if (strcmp(arg[i],"auto")==0)
schedule_user = 4;
if (strcmp(arg[i],"determine")==0)
schedule_user = 5;
if (schedule_user == 0)
error->all(FLERR,"Illegal pair_style command");
continue;
}
error->all(FLERR,"Illegal pair_style command");
}
if (nthreads < 0)
nthreads = comm->nthreads;
if (use_shared_arrays < 0) {
if (nthreads > 1 && atom->nlocal <= 2*nthreads)
use_shared_arrays = 1;
else use_shared_arrays = 0;
}
// check if running non-optimized code with
// optimization flags set
if (!use_optimized)
if (nthreads > 1 ||
use_shared_arrays ||
do_load_balance ||
schedule_user)
error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairSNAP::coeff(int narg, char **arg)
{
// read SNAP element names between 2 filenames
// nelements = # of SNAP elements
// elements = list of unique element names
if (narg < 6) error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
if (nelements) {
for (int i = 0; i < nelements; i++)
delete[] elements[i];
delete[] elements;
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(coeffelem);
}
nelements = narg - 4 - atom->ntypes;
if (nelements < 1) error->all(FLERR,"Incorrect args for pair coefficients");
char* type1 = arg[0];
char* type2 = arg[1];
char* coefffilename = arg[2];
char** elemlist = &arg[3];
char* paramfilename = arg[3+nelements];
char** elemtypes = &arg[4+nelements];
// insure I,J args are * *
if (strcmp(type1,"*") != 0 || strcmp(type2,"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
elements = new char*[nelements];
for (int i = 0; i < nelements; i++) {
char* elemname = elemlist[i];
int n = strlen(elemname) + 1;
elements[i] = new char[n];
strcpy(elements[i],elemname);
}
// read snapcoeff and snapparam files
read_files(coefffilename,paramfilename);
if (!quadraticflag)
ncoeff = ncoeffall - 1;
else {
// ncoeffall should be (ncoeff+2)*(ncoeff+1)/2
// so, ncoeff = floor(sqrt(2*ncoeffall))-1
ncoeff = sqrt(2*ncoeffall)-1;
ncoeffq = (ncoeff*(ncoeff+1))/2;
int ntmp = 1+ncoeff+ncoeffq;
if (ntmp != ncoeffall) {
printf("ncoeffall = %d ntmp = %d ncoeff = %d \n",ncoeffall,ntmp,ncoeff);
error->all(FLERR,"Incorrect SNAP coeff file");
}
}
// read args that map atom types to SNAP elements
// map[i] = which element the Ith atom type is, -1 if not mapped
// map[0] is not used
for (int i = 1; i <= atom->ntypes; i++) {
char* elemname = elemtypes[i-1];
int jelem;
for (jelem = 0; jelem < nelements; jelem++)
if (strcmp(elemname,elements[jelem]) == 0)
break;
if (jelem < nelements)
map[i] = jelem;
else if (strcmp(elemname,"NULL") == 0) map[i] = -1;
else error->all(FLERR,"Incorrect args for pair coefficients");
}
// clear setflag since coeff() called once with I,J = * *
int n = atom->ntypes;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
sna = new SNA*[nthreads];
// allocate memory for per OpenMP thread data which
// is wrapped into the sna class
#if defined(_OPENMP)
#pragma omp parallel default(none)
#endif
{
int tid = omp_get_thread_num();
sna[tid] = new SNA(lmp,rfac0,twojmax,
diagonalstyle,use_shared_arrays,
rmin0,switchflag,bzeroflag);
if (!use_shared_arrays)
sna[tid]->grow_rij(nmax);
}
printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff);
if (ncoeff != sna[0]->ncoeff) {
printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff);
error->all(FLERR,"Incorrect SNAP parameter file");
}
// Calculate maximum cutoff for all elements
rcutmax = 0.0;
for (int ielem = 0; ielem < nelements; ielem++)
rcutmax = MAX(2.0*radelem[ielem]*rcutfac,rcutmax);
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairSNAP::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style SNAP requires newton pair on");
// need a full neighbor list
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
#if defined(_OPENMP)
#pragma omp parallel default(none)
#endif
{
int tid = omp_get_thread_num();
sna[tid]->init();
}
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSNAP::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
return (radelem[map[i]] +
radelem[map[j]])*rcutfac;
}
/* ---------------------------------------------------------------------- */
void PairSNAP::read_files(char *coefffilename, char *paramfilename)
{
// open SNAP coefficient file on proc 0
FILE *fpcoeff;
if (comm->me == 0) {
fpcoeff = force->open_potential(coefffilename);
if (fpcoeff == NULL) {
char str[128];
sprintf(str,"Cannot open SNAP coefficient file %s",coefffilename);
error->one(FLERR,str);
}
}
char line[MAXLINE],*ptr;
int eof = 0;
int n;
int nwords = 0;
while (nwords == 0) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpcoeff);
if (ptr == NULL) {
eof = 1;
fclose(fpcoeff);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
}
if (nwords != 2)
error->all(FLERR,"Incorrect format in SNAP coefficient file");
// words = ptrs to all words in line
// strip single and double quotes from words
char* words[MAXWORD];
int iword = 0;
words[iword] = strtok(line,"' \t\n\r\f");
iword = 1;
words[iword] = strtok(NULL,"' \t\n\r\f");
int nelemfile = atoi(words[0]);
ncoeffall = atoi(words[1]);
// Set up element lists
memory->create(radelem,nelements,"pair:radelem");
memory->create(wjelem,nelements,"pair:wjelem");
memory->create(coeffelem,nelements,ncoeffall,"pair:coeffelem");
int *found = new int[nelements];
for (int ielem = 0; ielem < nelements; ielem++)
found[ielem] = 0;
// Loop over elements in the SNAP coefficient file
for (int ielemfile = 0; ielemfile < nelemfile; ielemfile++) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpcoeff);
if (ptr == NULL) {
eof = 1;
fclose(fpcoeff);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof)
error->all(FLERR,"Incorrect format in SNAP coefficient file");
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
nwords = atom->count_words(line);
if (nwords != 3)
error->all(FLERR,"Incorrect format in SNAP coefficient file");
iword = 0;
words[iword] = strtok(line,"' \t\n\r\f");
iword = 1;
words[iword] = strtok(NULL,"' \t\n\r\f");
iword = 2;
words[iword] = strtok(NULL,"' \t\n\r\f");
char* elemtmp = words[0];
double radtmp = atof(words[1]);
double wjtmp = atof(words[2]);
// skip if element name isn't in element list
int ielem;
for (ielem = 0; ielem < nelements; ielem++)
if (strcmp(elemtmp,elements[ielem]) == 0) break;
if (ielem == nelements) {
if (comm->me == 0)
for (int icoeff = 0; icoeff < ncoeffall; icoeff++)
ptr = fgets(line,MAXLINE,fpcoeff);
continue;
}
// skip if element already appeared
if (found[ielem]) {
if (comm->me == 0)
for (int icoeff = 0; icoeff < ncoeffall; icoeff++)
ptr = fgets(line,MAXLINE,fpcoeff);
continue;
}
found[ielem] = 1;
radelem[ielem] = radtmp;
wjelem[ielem] = wjtmp;
if (comm->me == 0) {
if (screen) fprintf(screen,"SNAP Element = %s, Radius %g, Weight %g \n",
elements[ielem], radelem[ielem], wjelem[ielem]);
if (logfile) fprintf(logfile,"SNAP Element = %s, Radius %g, Weight %g \n",
elements[ielem], radelem[ielem], wjelem[ielem]);
}
for (int icoeff = 0; icoeff < ncoeffall; icoeff++) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpcoeff);
if (ptr == NULL) {
eof = 1;
fclose(fpcoeff);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof)
error->all(FLERR,"Incorrect format in SNAP coefficient file");
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
nwords = atom->count_words(line);
if (nwords != 1)
error->all(FLERR,"Incorrect format in SNAP coefficient file");
iword = 0;
words[iword] = strtok(line,"' \t\n\r\f");
coeffelem[ielem][icoeff] = atof(words[0]);
}
}
// set flags for required keywords
rcutfacflag = 0;
twojmaxflag = 0;
// Set defaults for optional keywords
rfac0 = 0.99363;
rmin0 = 0.0;
diagonalstyle = 3;
switchflag = 1;
bzeroflag = 1;
quadraticflag = 0;
// open SNAP parameter file on proc 0
FILE *fpparam;
if (comm->me == 0) {
fpparam = force->open_potential(paramfilename);
if (fpparam == NULL) {
char str[128];
sprintf(str,"Cannot open SNAP parameter file %s",paramfilename);
error->one(FLERR,str);
}
}
eof = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpparam);
if (ptr == NULL) {
eof = 1;
fclose(fpparam);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
if (nwords != 2)
error->all(FLERR,"Incorrect format in SNAP parameter file");
// words = ptrs to all words in line
// strip single and double quotes from words
char* keywd = strtok(line,"' \t\n\r\f");
char* keyval = strtok(NULL,"' \t\n\r\f");
if (comm->me == 0) {
if (screen) fprintf(screen,"SNAP keyword %s %s \n",keywd,keyval);
if (logfile) fprintf(logfile,"SNAP keyword %s %s \n",keywd,keyval);
}
if (strcmp(keywd,"rcutfac") == 0) {
rcutfac = atof(keyval);
rcutfacflag = 1;
} else if (strcmp(keywd,"twojmax") == 0) {
twojmax = atoi(keyval);
twojmaxflag = 1;
} else if (strcmp(keywd,"rfac0") == 0)
rfac0 = atof(keyval);
else if (strcmp(keywd,"rmin0") == 0)
rmin0 = atof(keyval);
else if (strcmp(keywd,"diagonalstyle") == 0)
diagonalstyle = atoi(keyval);
else if (strcmp(keywd,"switchflag") == 0)
switchflag = atoi(keyval);
else if (strcmp(keywd,"bzeroflag") == 0)
bzeroflag = atoi(keyval);
else if (strcmp(keywd,"quadraticflag") == 0)
quadraticflag = atoi(keyval);
else
error->all(FLERR,"Incorrect SNAP parameter file");
}
if (rcutfacflag == 0 || twojmaxflag == 0)
error->all(FLERR,"Incorrect SNAP parameter file");
delete[] found;
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double PairSNAP::memory_usage()
{
double bytes = Pair::memory_usage();
int n = atom->ntypes+1;
bytes += n*n*sizeof(int);
bytes += n*n*sizeof(double);
bytes += 3*nmax*sizeof(double);
bytes += nmax*sizeof(int);
bytes += (2*ncoeffall)*sizeof(double);
bytes += (ncoeff*3)*sizeof(double);
bytes += sna[0]->memory_usage()*nthreads;
return bytes;
}

Event Timeline

c4science · Help