From 71ee2ecaa1e023ce358d6fa06f0f8c748dde7867 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 9 Nov 2016 14:52:39 -0500 Subject: [PATCH] integrate pair style tersoff/mod/c contributed by Ganga P Purja Pun (GMU) This includes docs, added testing and inclusion of USER-OMP support. --- doc/src/Eqs/pair_tersoff_mod_c.jpg | Bin 0 -> 4202 bytes doc/src/Eqs/pair_tersoff_mod_c.tex | 10 + doc/src/pair_tersoff_mod.txt | 38 ++-- examples/threebody/in.threebody | 32 +++ src/MANYBODY/pair_tersoff.h | 1 + src/MANYBODY/pair_tersoff_mod.h | 2 +- src/MANYBODY/pair_tersoff_mod_c.cpp | 196 ++++++++++++++++++ src/MANYBODY/pair_tersoff_mod_c.h | 66 +++++++ src/USER-OMP/pair_tersoff_mod_c_omp.cpp | 253 ++++++++++++++++++++++++ src/USER-OMP/pair_tersoff_mod_c_omp.h | 43 ++++ 10 files changed, 628 insertions(+), 13 deletions(-) create mode 100644 doc/src/Eqs/pair_tersoff_mod_c.jpg create mode 100644 doc/src/Eqs/pair_tersoff_mod_c.tex create mode 100644 src/MANYBODY/pair_tersoff_mod_c.cpp create mode 100644 src/MANYBODY/pair_tersoff_mod_c.h create mode 100644 src/USER-OMP/pair_tersoff_mod_c_omp.cpp create mode 100644 src/USER-OMP/pair_tersoff_mod_c_omp.h diff --git a/doc/src/Eqs/pair_tersoff_mod_c.jpg b/doc/src/Eqs/pair_tersoff_mod_c.jpg new file mode 100644 index 0000000000000000000000000000000000000000..2d2664eb5c880eb9ec358793765a4ef4fd8ebba7 GIT binary patch literal 4202 zcmbW4XH?V8wuk>AKHkRnBTjSzZ=ND)MOM?pb)hkyu#P<;WZ ziii|L@1aYTuEOCx=d63zx}Wa7d)A&Wv)AnXteN%OGZ*6*^8nKWHIy0v0s#Q%asU@o zfC@lGNeQ8(pn^ak)YMcobWnOaT3R{|7FI?mH|I4TZcZ+)>wI?wuHO)Xb8!hu+!DGc zfQ+ASb_!4!&Fm$eAfvctw;b zS@monH$0(;u#|i%xN=PgoBq%)pQx?Zb7~s)D=-dD{+j}Vw{9cF#3dx9q*d;#s;O(B z9y~HIG%|i{VrplPe&XPWaq{-@_45zFJ`H~n5g8R76Px<-Ra*M%HyN1)g+;|BrDf$6 zwRQM(CYATs|;2Q!mX@QP5fDCt3LJXvob!l{*|9_GF7ubJuO#*aa(B`|zdKX^g6$Z<9bZZQaeND>CtR z=UV%T zL`phDi2LmQRk`M+#I%=x%55$HRaV0bpp8BEcY-4K4BN|v^@YZfnFXSo7@CE7^=l38 zR{N<(bU62!Lr#5_A_dgYRpm$9#MsRvNmBjeiDDEKkp>pHzWY=cB2Z$nDbGjrrESr^ zZ)$=$>!@Rf9aG!HD86o+#ZYrdwOayG0`>107v5!;I<_&oh`o(7(=^!@&e>mC?KkA* zv-^i|{1%NmPpHv$cKOESa)A@OuGz%+XNjpe{F8M1AKgkkkLA?-xtNm~p`dziZ?^2B zF=GZ6m|juJpsFrHhjiy3oysdqn%ermJX%Jf4~JM&Xn6GRDNH;X#Ieh5zTH#9`Bwz- z?@=s&FCCvLXbkOa^E5CPsS$`LQ7S0k!MEan8OBZGPBl)9ZYzYty^^%zb6~LJBrDhXmO+9XAKVyVHq@KZvI9p9zW#rg?;Wf+xO-+SMpOizaqzOH3_A zO6ui}_m~7tBP6b>O?ziGNWpr8#^jcC@BhW=dwB<&r>{O^^O!qXcLFw?Gk56|s;Zyk zWAeh_=_{Y99ela^mMI^fWv)BZb|~!M2zYkhLA(IYgdNFBqYq^|7quBXk3_LIPjuwE zaL4I_$<|NfTk`rk+9lytuIuc_giO^+&a}l@%U&3(`E3`5*~Tq%G$F3oY3_<|25Cfy zq^+Ry{N@~4_Q>ZgUhiJv`bRmfzs8gn>D`XCzaMBz_guxIM^^)CZrLnXz-1J;729El zHj3O7)@IWv+)2S6IK{3gEzYTdWx6&>=7Tc)CPb@7&cMb+Y@M{HJs4Fjy!7k7c#NMw zgKeK));l73P*qAY!b^1BP?pBuWYevaVRGhZ=c9h%1;9CdE&e8_VPDt&KFiGwoNd?(aW{RF!}crEia z!+or6dknIVJ3cn;k^s5*v;>yP%kFSD<7*j1gBYrtvHUiitFnT;@sIJHw{M|kVWWA2SD4vrgUe;^diB>5PY-TpFVL2rPwyq;RXHu{&r~q%e%I4%^=?+M z2Zqn&eMg8Gn}oYs)^CF5UVnja@AOu4^e*Tt4u5EMm38JSR`M~LW~=FgfMt){^B5jk zw2dAzvX@34wA46$D((;F8y8p=g1_Lh-uilUojE(tr%z71$h#X&?k^g%%_3%krfv9Q zE#>0(X_f8ciIwVbbx`o)fG#E-E&MJI-{8(0f`Z4(#f=4ZaU$d zZ#U;{C5ok+f0TDtO!kh?b%$xr!4nV7Aju&vYmK9hDR^P3W=PRDh@mQWHmtCtEMTfD zrEf<#;gvX1tF1$>r^i37SIu(;^fE{1eXkW7{U^L={me3_$}9Ot^wqaXx&4q8FAskv z8PD}@z-h>E`$?vFgzEyzf&}_rJW;O@69thx*x9(<4XF5g5&Ry_Nb>V06j^`J0Lz9L zm`)4NWL*)YGAT3x7#;t|`|Uh!2x>>Ir}`t2`Dm>TP>WLzQG}Mrip# z7@aQ*LOjuP`5WYsD%-~ZRmKFoSuXR&gX*JGehc|qy3&)1m3kURuvrZu)Mm0YAP8=$ z;QWwv+#Fa{)I-~M0bt$l2ADv!JDYP}_>&X0Ber^Sa~a2ZN4JqU$%L+B(NlCk81GRk zq|EXZl(N~=lTC9v16XVM~lULe)9(!H@&lpd}zRN=>e1a`S&dbHJmF&1fe@-VWusD-8fL{=+ z-TJ7QyZDUi&U+(nA;ZzqQj==LBNCT?kF?<*L-4wUQq?Spsvk>*ea-dbw`aWnfO$cv z6eELKFM6ITvI<&}-Zy4!>GXPuL zqudr26asBNPrDTPbY_d5eN|}`dCh4h+tNa88ibssP`)tnMeLd)08j##4fgLQ`yX!{ zyy|R)8997LJKn8uRQ0`}lls;RcgZ*BkEY9}d%1zH&UKeHSzQV6a=ri}kN>(h++)uk zLwBpiZTD@6*3tA8_a#%k&Ex?`oE^q%jNc=D&}|PY7f!MP4v&^(9c*KP!S6nJq7UuB z?X57g7Ki6*hP|`;)Zhm5?{YMYuEw33t^r#N`OqiBKW$D6-_VJ094@>%rel0MT_Pk> zIf&>fPzOZA`axOtzw9qX=aK5LD_)qTABCSe3Y0X03Xev>J{`+xd!M0;YJyb0yV0k8 zyw9s)n_ND=P$!hiE-Ufn6O5yeDv)Lqr(n^LxT(?_X)#52Kl`1*ULT1mP?wS+&XqMt ziCENN@qXM?$?>6=eOkZPa*}ss(oZpW-K1Z3-Dj?~3B)s$^QUnh@-GoLbl#2m!wIaT zMu|~&?&uPn5)ZLS!7(F_9m8i$ZL7m~a0USLG;LfW$}m;Vl9jttYKY!)`2*mLgp_ld zQa zzrABo1qJn-K>MA}1`-B}yB^&w^NP@m$AyBi)WO!|(wvr4PaZ zoN+Ps>5?0Xd#q0e*C#Y@J)h){6j`ePuI1sNdWsyAA=&&Tr8o;2obZqJ+o12tT`_}; zC3wezP((En+4!BI*xg?)dur-WCW@EFUk&ztA&i5Z!1znO%KQBG?$(Ef=afk$CJ$_} z>|}OEDTO9xE+fqSqOJ~7B+}|^i>V+ues2`zKR+ncJ~T5Ieef%QBqdVcDP(`E~F^ z7tM1Qa!Y>d=C`_;;xIw9*XJTXUcdN2gOpF=Ey4UG zeai_Y>j5rJ;kR&&?hp)FUlSbS*nWOo<>R9Z;ocn*!0+NP4eZAr_1Vjxd7B6US}7UX zDWTQ~nv{`Z^W9hK8k&QQnjLjHDA&5*e#25e@IS*}qW%neS_q5eCoH2(#Sqtd%@YkS z#XX4ebExQn1;dnz=7QV4jODCmps0ot|MbB;0{}2u-(fYwr|WP#dRL`&r6}|GXF!wf zdtXJH6_lBvy|Wy+2ieAS)s!y30eu^UA)P%&rG}ViSXi!pW6;ltc*VG!PEnCwomaW0 zo4B0acWDj|F2_t=vmx>WD=4@ORxg#=t1F8xx1294bIp;AQ0Ow2iJqyrXPkK@0kuZ0 zuTUi%(PD}r627V+?-ki6z6hmwTo6jrU9(GiiY~9MX=`bCKf6s>ucK~?L*QJeukbBp z#qTNt^c8~Ogqj8&x_~;X#?)CG2duj*f5S^!RAWzk0<2v@g2r=ndr4AEBJ#$|jyZbd zcs-wjBOY&sfG)>Z^(6#9uA_i$I9HaS==AZfd`>Ikvk{Sg1ggFRNr9hvbV>9WueEbn z1HQ3$GPZS$k1HG*Z>hzz +#include +#include +#include +#include "pair_tersoff_mod_c.h" +#include "atom.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "force.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +#include "math_const.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define MAXLINE 1024 +#define DELTA 4 + +/* ---------------------------------------------------------------------- */ + +void PairTersoffMODC::read_file(char *file) +{ + int params_per_line = 21; + char **words = new char*[params_per_line+1]; + + memory->sfree(params); + params = NULL; + nparams = maxparam = 0; + + // open file on proc 0 + + FILE *fp; + if (comm->me == 0) { + fp = force->open_potential(file); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open Tersoff potential file %s",file); + error->one(FLERR,str); + } + } + + // read each line out of file, skipping blank lines or leading '#' + // store line of params if all 3 element tags are in element list + + int n,nwords,ielement,jelement,kelement; + char line[MAXLINE],*ptr; + int eof = 0; + + while (1) { + if (comm->me == 0) { + ptr = fgets(line,MAXLINE,fp); + if (ptr == NULL) { + eof = 1; + fclose(fp); + } 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; + + // concatenate additional lines until have params_per_line words + + while (nwords < params_per_line) { + n = strlen(line); + if (comm->me == 0) { + ptr = fgets(&line[n],MAXLINE-n,fp); + if (ptr == NULL) { + eof = 1; + fclose(fp); + } 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); + if ((ptr = strchr(line,'#'))) *ptr = '\0'; + nwords = atom->count_words(line); + } + + if (nwords != params_per_line) + error->all(FLERR,"Incorrect format in Tersoff potential file"); + + // words = ptrs to all words in line + + nwords = 0; + words[nwords++] = strtok(line," \t\n\r\f"); + while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue; + + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next line + + for (ielement = 0; ielement < nelements; ielement++) + if (strcmp(words[0],elements[ielement]) == 0) break; + if (ielement == nelements) continue; + for (jelement = 0; jelement < nelements; jelement++) + if (strcmp(words[1],elements[jelement]) == 0) break; + if (jelement == nelements) continue; + for (kelement = 0; kelement < nelements; kelement++) + if (strcmp(words[2],elements[kelement]) == 0) break; + if (kelement == nelements) continue; + + // load up parameter settings and error check their values + + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); + } + + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].powerm = atof(words[3]); + params[nparams].lam3 = atof(words[4]); + params[nparams].h = atof(words[5]); + params[nparams].powern = atof(words[6]); + params[nparams].beta = atof(words[7]); + params[nparams].lam2 = atof(words[8]); + params[nparams].bigb = atof(words[9]); + params[nparams].bigr = atof(words[10]); + params[nparams].bigd = atof(words[11]); + params[nparams].lam1 = atof(words[12]); + params[nparams].biga = atof(words[13]); + params[nparams].powern_del = atof(words[14]); + params[nparams].c1 = atof(words[15]); + params[nparams].c2 = atof(words[16]); + params[nparams].c3 = atof(words[17]); + params[nparams].c4 = atof(words[18]); + params[nparams].c5 = atof(words[19]); + params[nparams].c0 = atof(words[20]); + + // currently only allow m exponent of 1 + + params[nparams].powermint = int(params[nparams].powerm); + + if ( + params[nparams].lam3 < 0.0 || params[nparams].powern < 0.0 || + params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 || + params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 || + params[nparams].bigd < 0.0 || + params[nparams].bigd > params[nparams].bigr || + params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0 || + params[nparams].powerm - params[nparams].powermint != 0.0 || + (params[nparams].powermint != 3 && params[nparams].powermint != 1)) + error->all(FLERR,"Illegal Tersoff parameter"); + + nparams++; + } + + delete [] words; +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoffMODC::repulsive(Param *param, double rsq, double &fforce, + int eflag, double &eng) +{ + double r,tmp_fc,tmp_fc_d,tmp_exp; + + r = sqrt(rsq); + tmp_fc = ters_fc(r,param); + tmp_fc_d = ters_fc_d(r,param); + tmp_exp = exp(-param->lam1 * r); + fforce = -param->biga * tmp_exp * (tmp_fc_d - tmp_fc * param->lam1) / r - param->c0 * tmp_fc_d / r; + if (eflag) eng = tmp_fc * param->biga * tmp_exp + param->c0 * tmp_fc; +} + diff --git a/src/MANYBODY/pair_tersoff_mod_c.h b/src/MANYBODY/pair_tersoff_mod_c.h new file mode 100644 index 0000000000..5372bc2b15 --- /dev/null +++ b/src/MANYBODY/pair_tersoff_mod_c.h @@ -0,0 +1,66 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(tersoff/mod/c,PairTersoffMODC) + +#else + +#ifndef LMP_PAIR_TERSOFF_MOD_C_H +#define LMP_PAIR_TERSOFF_MOD_C_H + +#include "pair_tersoff_mod.h" + +namespace LAMMPS_NS { + +class PairTersoffMODC : public PairTersoffMOD { + public: + PairTersoffMODC(class LAMMPS *lmp) : PairTersoffMOD(lmp) {}; + ~PairTersoffMODC() {} + + protected: + void read_file(char *); + void repulsive(Param *, double, double &, int, double &); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Cannot open Tersoff potential file %s + +The specified potential file cannot be opened. Check that the path +and name are correct. + +E: Incorrect format in Tersoff potential file + +Incorrect number of words per line in the potential file. + +E: Illegal Tersoff parameter + +One or more of the coefficients defined in the potential file is +invalid. + +E: Potential file has duplicate entry + +The potential file has more than one entry for the same element. + +E: Potential file is missing an entry + +The potential file does not have a needed entry. + +*/ diff --git a/src/USER-OMP/pair_tersoff_mod_c_omp.cpp b/src/USER-OMP/pair_tersoff_mod_c_omp.cpp new file mode 100644 index 0000000000..340eb3ebc5 --- /dev/null +++ b/src/USER-OMP/pair_tersoff_mod_c_omp.cpp @@ -0,0 +1,253 @@ +/* ---------------------------------------------------------------------- + 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 +#include "pair_tersoff_mod_c_omp.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" + +#include "suffix.h" +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairTersoffMODCOMP::PairTersoffMODCOMP(LAMMPS *lmp) : + PairTersoffMODC(lmp), ThrOMP(lmp, THR_PAIR) +{ + suffix_flag |= Suffix::OMP; + respa_enable = 0; +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoffMODCOMP::compute(int eflag, int vflag) +{ + if (eflag || vflag) { + ev_setup(eflag,vflag); + } else evflag = vflag_fdotr = vflag_atom = 0; + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = list->inum; + +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(eflag,vflag) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr); + + if (evflag) { + if (eflag) { + if (vflag_atom) eval<1,1,1>(ifrom, ito, thr); + else eval<1,1,0>(ifrom, ito, thr); + } else { + if (vflag_atom) eval<1,0,1>(ifrom, ito, thr); + else eval<1,0,0>(ifrom, ito, thr); + } + } else eval<0,0,0>(ifrom, ito, thr); + + thr->timer(Timer::PAIR); + reduce_thr(this, eflag, vflag, thr); + } // end of omp parallel region +} + +template +void PairTersoffMODCOMP::eval(int iifrom, int iito, ThrData * const thr) +{ + int i,j,k,ii,jj,kk,jnum; + tagint itag,jtag; + int itype,jtype,ktype,iparam_ij,iparam_ijk; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,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; + + evdwl = 0.0; + + const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; + dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; + const tagint * _noalias const tag = atom->tag; + const int * _noalias const type = atom->type; + const int nlocal = atom->nlocal; + + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + double fxtmp,fytmp,fztmp; + + // 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].x; + ytmp = x[i].y; + ztmp = x[i].z; + fxtmp = fytmp = fztmp = 0.0; + + // two-body interactions, skip half of them + + 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].z < ztmp) continue; + if (x[j].z == ztmp && x[j].y < ytmp) continue; + if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue; + } + + jtype = map[type[j]]; + + delx = xtmp - x[j].x; + dely = ytmp - x[j].y; + delz = ztmp - x[j].z; + rsq = delx*delx + dely*dely + delz*delz; + + iparam_ij = elem2param[itype][jtype][jtype]; + if (rsq > params[iparam_ij].cutsq) continue; + + repulsive(¶ms[iparam_ij],rsq,fpair,EFLAG,evdwl); + + fxtmp += delx*fpair; + fytmp += dely*fpair; + fztmp += delz*fpair; + f[j].x -= delx*fpair; + f[j].y -= dely*fpair; + f[j].z -= delz*fpair; + + if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1, + evdwl,0.0,fpair,delx,dely,delz,thr); + } + + // three-body interactions + // skip immediately if I-J is not within cutoff + double fjxtmp,fjytmp,fjztmp; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtype = map[type[j]]; + iparam_ij = elem2param[itype][jtype][jtype]; + + delr1[0] = x[j].x - xtmp; + delr1[1] = x[j].y - ytmp; + delr1[2] = x[j].z - ztmp; + rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + 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; + + 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].x - xtmp; + delr2[1] = x[k].y - ytmp; + delr2[2] = x[k].z - ztmp; + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + if (rsq2 > params[iparam_ijk].cutsq) continue; + + zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2); + } + + // pairwise force due to zeta + + force_zeta(¶ms[iparam_ij],rsq1,zeta_ij,fpair,prefactor,EFLAG,evdwl); + + 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],thr); + + // attractive term via loop over k + + 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].x - xtmp; + delr2[1] = x[k].y - ytmp; + delr2[2] = x[k].z - ztmp; + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + if (rsq2 > params[iparam_ijk].cutsq) continue; + + attractive(¶ms[iparam_ijk],prefactor, + rsq1,rsq2,delr1,delr2,fi,fj,fk); + + fxtmp += fi[0]; + fytmp += fi[1]; + fztmp += fi[2]; + fjxtmp += fj[0]; + fjytmp += fj[1]; + fjztmp += fj[2]; + f[k].x += fk[0]; + f[k].y += fk[1]; + f[k].z += fk[2]; + + if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr); + } + f[j].x += fjxtmp; + f[j].y += fjytmp; + f[j].z += fjztmp; + } + f[i].x += fxtmp; + f[i].y += fytmp; + f[i].z += fztmp; + } +} + +/* ---------------------------------------------------------------------- */ + +double PairTersoffMODCOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairTersoffMOD::memory_usage(); + + return bytes; +} diff --git a/src/USER-OMP/pair_tersoff_mod_c_omp.h b/src/USER-OMP/pair_tersoff_mod_c_omp.h new file mode 100644 index 0000000000..b3891364ca --- /dev/null +++ b/src/USER-OMP/pair_tersoff_mod_c_omp.h @@ -0,0 +1,43 @@ +/* -*- 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(tersoff/mod/c/omp,PairTersoffMODCOMP) + +#else + +#ifndef LMP_PAIR_TERSOFF_MOD_C_OMP_H +#define LMP_PAIR_TERSOFF_MOD_C_OMP_H + +#include "pair_tersoff_mod_c.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairTersoffMODCOMP : public PairTersoffMODC, public ThrOMP { + + public: + PairTersoffMODCOMP(class LAMMPS *); + + virtual void compute(int, int); + virtual double memory_usage(); + + private: + template + void eval(int ifrom, int ito, ThrData * const thr); +}; + +} + +#endif +#endif