From 71ee2ecaa1e023ce358d6fa06f0f8c748dde7867 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 9 Nov 2016 14:52:39 -0500 Subject: [PATCH 1/3] 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 From c55fd502e0b9bda547d97e48056aa4333f53b955 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 9 Nov 2016 15:04:24 -0500 Subject: [PATCH 2/3] correct typo in formula --- doc/src/Eqs/pair_tersoff_mod_c.jpg | Bin 4202 -> 4211 bytes doc/src/Eqs/pair_tersoff_mod_c.tex | 2 +- doc/src/pair_tersoff_mod.txt | 1 + 3 files changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/src/Eqs/pair_tersoff_mod_c.jpg b/doc/src/Eqs/pair_tersoff_mod_c.jpg index 2d2664eb5c880eb9ec358793765a4ef4fd8ebba7..311ccc81ebaa809dea7da497439a7c79a405acd3 100644 GIT binary patch delta 1980 zcmV;t2SfPkAoC!w&Io_kKdxOqiKtlWq1)_pPaU1&6Y_xVwmAo(4pmR8w?A}N#mv)2 z(5<{FE4+X@w7WK^@i5@$9-aL~bKepH{{V$|tup=fgL;3y7#00Mpbn;rDRxxm-tuue?IK zMB*(D(bxM%{E~G7yqmdlKh^d*8SjkoU6r%XCA4pI3d1oOc+vNTK7p_^`A`J3R~KP3 z_^-yg+*gacNYP$9SZ$k>nsU3H+;`dvg$E<4>x^BtKk$v+#iQFfy}Q(|C%T3?2!=Q! zUAY5loGR`>02zNJi9O8#ZfaV6rl(}~_qOrR3Hi22GBP?E_}lwJ_x^tFKkwoH02=!P zN!Alsw7eoeDtk%fwYZCR;F2$);~;{ha5y7@UmMT|*gF3Jg{(YDb)?OH8HSs4_LM=< z&eArww;=mv9LofHAe_`AcCU;H<{@cx!FAhbwfP1|>b z+wy0&;Jbes(vMt*YeicdOP0ASgtmKkh@tEsBM1ktw;AI*6i^0!tEHy1pk8TLYzp1m zJQ6alKJt)A9lGYJ=vNbJz8YI=$7pp3qc)eD56yoZz+s5buMDTSV0+f(F=35nd6qR2 zI02jvSPbJQ+b8p&3YspHZJ}tET1}?!v)RJ1p>Cuo#!tEGYoF6CC5z%7t8i?t5+Mph z6opR54Vh$+vEjaENWjiOL)cecS^Dj~cUm3fGVa)7B<*sijz=Ja+JG>0TWFi&4z(q^ z?S_BC(hEX!lO(QMXDoZ1G=G_^dW3)59xu~vly0}uOI}XE@Frv_5`A{NxgD^cqO^27 z8+$D^S_t;F$@@giShGnZ080{iC#TQ=2LM(0;ffH#DAhMHcG@wPJxTYV2iGB?>sHzq z!=DXY-A6OWC9R|%x=q_lqK(biDC$Y+U44I1({z0|M~ULoY@oG`&RJqjpwI>%g6$u86*Z$p6akfap{2&Cbt`+h?C##*?9xdHQ@Sw4hYSeFC9r?mJ!_xwSBGuwybv^f zGRjjE$gmR{w2d0d+h&-KS(|G)&nG87_1e)u9A}QAOTPwrW(JP#{70wT$05m$t{Gz| z({iaNJ+OJ?=X+MyY_t!srQm;_q}VQqJTK>hKW6=gq>1mA&Zz@XB zsgVH+%0a%*dq_TMHL!E4{ zB#IF8yvU@3COo7<-4%aA%PXAYIPQ3^x`HwY^`H(%Rkn(Mh z*}TalXKyMpLM~X7fs^$llapJ^3`r}>v8a+j1GpTp9Ws5ge>#5veH2kZ8QLbDCYNC? zysoa#d3A6K^u)8tARka11MVwEiYNmibO_+ouWxK`E-c>pmF^ZWge=!18L_TT{}m= z@et9h0_av{KF=5obrLc*j&bGC2UA5A!gxgwjkH^z5$Z~>_W6G$T|OyLI@&<2QAk!8 z-bgG->ySDCC*0C{gNcv%Nc!@8b?cLdf=h3&?lpAXm})b2E^ zIE}>jGx=r91GI3bEJ6}`0g>y*#*eDoYnnBlsb<96UBa_RxL|<8I6mW{pb8dq*m%Ff zw%1oSi16G-9pHbIZt~O@ln#IAGe~xv5TV(C2dG~|S9}R= zHLrp+p>Q_OG)Cq{{_9G{Lcf80eiXX?m1n7HmJwaX@yXW8t7B<$*ClY4 z&u;M)J%i+70QL6cJZEBx0MFHQwAOSB%_{AITf1k1Mpeh&QV9dMT-7}a;%!&MOKo`V z4xtog((_^YqnHda8TH|m_Y4nu+@>rstj{vWqDOxKGl9zioMiiC{&WFTMbd4wEfUK~ zwB7c5I93!b)P)$y_dQK>`emfCd_&c44VA(~AxL77so3GOERq&HH_XWx8OR8G3hS#s zUAKV$xw{1& zNj)pCsyc3urs(lJT5XgTvAN4ENwgXO!|;Dyq@}gJZgnBGZ9(RXJxmCqIFI6~xO*_F zlR$O9iOAF2Y6)eg>CYT9c~e4-JdzFjWB>z&jeFqk?_I4F0mk?g%WC?~jm@RRQ!b}- z<-DZHLv&eaTgYuF|Aa&xd=+PY;#QMIUClG2f{{U!OG3Z#TtjYBZq5hJ48tH#& zl0zJ&43Vpt+j*O zCCvIamuGQ0i^O@J6f99A{peATx?V`*UEBHykWfz&&K3;hsf&$60>MbfO})^DzaY}3Is zi4(~qKb9mR-)>39Thp)upaQqdF&UO)8kto}tC7eD836hm56DmiRyJ_y+I_W+y6uA6 zC5AERssJCLs}xZ{8QOlGC8npVT|28S?CB~2)k4IH&*G-EXrh2JuQW8c)ee89Z+AVN z+uQw`Ng)b%Mi{W+ff)q05Bo>0bN(vu?Y*~x29Ks$N@95y0%KN@qgh*Q(-F%vZD%>; zWo}Nbse%I}^ad^{P?nR&#%Ny5@OfX&xe%F^7?vR#Ki`;H2Okqa8TO$0X;Xpblfe zdVR&0hZ^ScCW`Y?x|>j%HW(;na>}H3A$E*=RFTOxyL~0(QeN3XV9gv$Et?}q#eQi5 z$hjj2Xu$sWt#ouy2OaRyC$#XZX}1>laJ8M~$|1S7cubL*osocGG75i@)1_SSlt)NQ z&ksw%vz|$9v;20_Idmia(#(hb&v9KXF`#Lsdu9atJZ1)hv z-hMvK=1C(vc~O!Pa>SerpQ$9An&@Z(wxev9YyBi!!Ti+!uc&#*B>^z&|$%an$bPI5jthB9p=24K&!bBtK}hv|ln5 z<&3S!jeS_N1sLjUtfol}QnYO(YUVcHg@$&v2*~s}C-a~Tj}2*Z=>8bDw7AFGioDk( z9D_8BUQ8b2b2Te`dtYha>Gpb^skXkhj^YAAF~~A9I#xIQAYKyx0HJeC{{X;s{{RXA F|Jmii%;f+8 diff --git a/doc/src/Eqs/pair_tersoff_mod_c.tex b/doc/src/Eqs/pair_tersoff_mod_c.tex index 2bfa4ce44a..8cea2d382c 100644 --- a/doc/src/Eqs/pair_tersoff_mod_c.tex +++ b/doc/src/Eqs/pair_tersoff_mod_c.tex @@ -4,7 +4,7 @@ \begin{document} \begin{eqnarray*} - V_{ij} & = & f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij} + c_0) \right] + V_{ij} & = & f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) + c_0 \right] \end{eqnarray*} \end{document} diff --git a/doc/src/pair_tersoff_mod.txt b/doc/src/pair_tersoff_mod.txt index ad27af63d2..f74f322b1e 100644 --- a/doc/src/pair_tersoff_mod.txt +++ b/doc/src/pair_tersoff_mod.txt @@ -42,6 +42,7 @@ The summations in the formula are over all neighbors J and K of atom I within a cutoff distance = R + D. The {tersoff/mod/c} style differs from {tersoff/mod} only in the formulation of the V_ij term, where it contains an additional c0 term. + :c,image(Eqs/pair_tersoff_mod_c.jpg) From 3d3a99c082f723d18bebab0eb8d7d08817897b4c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 9 Nov 2016 16:50:34 -0500 Subject: [PATCH 3/3] added missing potential for tersoff/mod/c --- potentials/Si.tersoff.modc | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 potentials/Si.tersoff.modc diff --git a/potentials/Si.tersoff.modc b/potentials/Si.tersoff.modc new file mode 100644 index 0000000000..6c8224667f --- /dev/null +++ b/potentials/Si.tersoff.modc @@ -0,0 +1,10 @@ +# DATE: 2016-11-09 CONTRIBUTOR: Ganga P Purja Pun (George Mason University, Fairfax) CITATION: Unknown +# +# Format: +# element1 element2 element3 +# beta alpha h eta +# beta_ters lam2 B R D lam1 A +# n c1 c2 c3 c4 c5 C +Si Si Si 3.00000000 1.80536502 -0.38136087 2.16152496 +1 1.39343356 117.78072440 2.87478837 0.33090566 3.18011795 3198.51383127 +1.98633876 0.20123243 614230.04310619 996439.09714140 3.33560562 25.20963770 -0.00592042