diff --git a/src/Makefile b/src/Makefile
index 2d6cd12180..25184d998f 100755
--- a/src/Makefile
+++ b/src/Makefile
@@ -50,7 +50,7 @@ PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars user-cuda \
user-diffraction user-drude user-eff user-fep user-h5md \
user-intel user-lb \
user-misc user-molfile user-omp user-phonon user-qmmm user-qtb \
- user-quip user-reaxc user-smd user-sph user-tally
+ user-quip user-reaxc user-smd user-smtbq user-sph user-tally
PACKLIB = compress gpu kim kokkos meam poems python reax voronoi \
user-atc user-awpmd user-colvars user-cuda user-h5md user-molfile \
diff --git a/src/USER-SMTBQ/README b/src/USER-SMTBQ/README
new file mode 100644
index 0000000000..19fac74388
--- /dev/null
+++ b/src/USER-SMTBQ/README
@@ -0,0 +1,13 @@
+This package implements the Second Moment Tight Binding - QEq (SMTB-Q)
+potential for the description of the ionocovalent bond in oxides.
+
+Authors: Nicolas Salles, Emile Maras, Olivier Politano, Robert T\'e9tot
+
+Contact email: lammps@u-bourgogne.fr
+
+See the doc page for the pair_style smtbq command to get started.
+
+There are potential files for this potential in the potentials dir.
+
+There are example scripts for using this package in
+examples/USER/smtbq.
diff --git a/src/USER-SMTBQ/pair_smtbq.cpp b/src/USER-SMTBQ/pair_smtbq.cpp
new file mode 100644
index 0000000000..52b431a8c2
--- /dev/null
+++ b/src/USER-SMTBQ/pair_smtbq.cpp
@@ -0,0 +1,4241 @@
+/* ----------------------------------------------------------------------
+ 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.
+ ------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ The SMTBQ code has been developed with the financial support of CNRS and
+ of the Regional Council of Burgundy (Convention n¡ 2010-9201AAO037S03129)
+
+ Copyright (2015)
+ Universite de Bourgogne : Nicolas SALLES, Olivier POLITANO
+ Universite de Paris-Sud Orsay : R. Tetot
+ Aalto University (Finland) : E. Maras
+
+ Please cite the related publication:
+ N. Salles, O. Politano, E. Amzallag and R. Tetot,
+ Comput. Mater. Sci., 111 (2016) 181-189
+
+ Contact : lammps@u-bourgogne.fr
+
+ This program is free software; you can redistribute it and/or
+ modify it under the terms of the GNU General Public License as
+ published by the Free Software Foundation; either version 2 of
+ the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ See the GNU General Public License for more details:
+ .
+ ------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_smtbq.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "group.h"
+#include "update.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+#include "domain.h"
+
+#include
+#include
+#include
+#include
+
+//to DELETE
+//#include
+//#include
+//#include
+//#include
+//to DELETE
+using namespace std;
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define MAXLINE 2048
+#define MAXTOKENS 2048
+#define DELTA 4
+#define PGDELTA 1
+#define MAXNEIGH 24
+#define Pi 4*atan(1)
+
+/* ---------------------------------------------------------------------- */
+
+PairSMTBQ::PairSMTBQ(LAMMPS *lmp) : Pair(lmp)
+{
+ MPI_Comm_rank(world,&me);
+ MPI_Comm_size(world,&nproc);
+
+ single_enable = 0;
+ restartinfo = 0;
+ one_coeff = 1;
+
+ nmax = 0;
+ rmin = 0.0;
+ dr = 0.0;
+ ds = 0.0;
+ kmax = 0;
+
+ nelements = 0;
+ elements = NULL;
+ nparams = 0;
+ maxparam = 0;
+ params = NULL;
+ intparams = NULL;
+
+ intype = NULL;
+ coultype = NULL;
+ fafb = NULL;
+ dfafb = NULL;
+ potqn = NULL;
+ dpotqn = NULL;
+ Vself = 0.0;
+ tabsmb = NULL;
+ tabsmr = NULL;
+ dtabsmb = NULL;
+ dtabsmr = NULL;
+
+ sbcov = NULL;
+ coord = NULL;
+ sbmet = NULL;
+ chimet = NULL;
+ ecov = NULL;
+
+ potmad = NULL;
+ potself = NULL;
+ potcov = NULL;
+ qf = NULL;
+ q1 = NULL;
+ q2 = NULL;
+ tab_comm = NULL;
+
+ nvsm = NULL;
+ vsm = NULL;
+ flag_QEq = NULL;
+ nQEqaall = NULL;
+ nQEqcall = NULL;
+ nQEqall = NULL;
+ nteam = 0;
+ cluster = 0;
+
+ Nevery = 0.0;
+ Neverypot = 0.0;
+
+ fct = NULL;
+
+
+ maxpage = 0;
+
+ // set comm size needed by this Pair
+
+ comm_forward = 1;
+ comm_reverse = 1;
+}
+
+/* ----------------------------------------------------------------------
+ check if allocated, since class can be destructed when incomplete
+ ------------------------------------------------------------------------- */
+
+PairSMTBQ::~PairSMTBQ()
+{
+ int i;
+ if (elements) {
+ for ( i = 0; i < nelements; i++) delete [] elements[i];
+
+ for( i = 0; i < atom->ntypes ; i++ ) free( params[i].nom );
+ for( i = 1; i <= maxintparam ; i++ ) free( intparams[i].typepot );
+ for( i = 1; i <= maxintparam ; i++ ) free( intparams[i].mode );
+ }
+
+ free(QEqMode);
+ free(QInitMode);
+ free(writepot);
+ free(writeenerg);
+ free(Bavard);
+
+ delete [] elements;
+ memory->sfree(params);
+ memory->sfree(intparams);
+
+ memory->destroy(intype);
+ memory->destroy(coultype);
+ memory->destroy(fafb);
+ memory->destroy(dfafb);
+ memory->destroy(potqn);
+ memory->destroy(dpotqn);
+
+ memory->destroy(ecov);
+ memory->destroy(sbcov);
+ memory->destroy(coord);
+ memory->destroy(sbmet);
+ memory->destroy(tabsmb);
+ memory->destroy(tabsmr);
+ memory->destroy(dtabsmb);
+ memory->destroy(dtabsmr);
+
+ memory->destroy(potmad);
+ memory->destroy(potself);
+ memory->destroy(potcov);
+ memory->destroy(chimet);
+
+ memory->destroy(nvsm);
+ memory->destroy(vsm);;
+ memory->destroy(flag_QEq);
+
+ memory->destroy(nQEqall);
+ memory->destroy(nQEqcall);
+ memory->destroy(nQEqaall);
+
+ memory->destroy(qf);
+ memory->destroy(q1);
+ memory->destroy(q2);
+ memory->destroy(tab_comm);
+
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+ delete [] map;
+ delete [] esm;
+ }
+
+ memory->destroy(fct);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSMTBQ::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");
+
+ map = new int[n+1];
+ esm = new double[n];
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+ ------------------------------------------------------------------------- */
+
+void PairSMTBQ::settings(int narg, char **arg)
+{
+ if (narg > 0) error->all(FLERR,"Illegal pair_style command");
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+ ------------------------------------------------------------------------- */
+
+void PairSMTBQ::coeff(int narg, char **arg)
+{
+ int i,j,n;
+
+ if (!allocated) allocate();
+
+ if (narg != 3 + atom->ntypes)
+ error->all(FLERR,"Incorrect args for pair coefficients");
+
+ // insure I,J args are * *
+
+ if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
+ error->all(FLERR,"Incorrect args for pair coefficients");
+
+ // read args that map atom types to elements in potential file
+ // map[i] = which element the Ith atom type is, -1 if NULL
+ // nelements = # of unique elements
+ // elements = list of element names
+
+ if (elements) {
+ for (i = 0; i < nelements; i++) delete [] elements[i];
+ delete [] elements;
+ }
+ elements = new char*[atom->ntypes];
+ for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;
+
+ nelements = 0;
+ for (i = 3; i < narg; i++) {
+ if (strcmp(arg[i],"NULL") == 0) {
+ map[i-2] = -1;
+ continue;
+ }
+ for (j = 0; j < nelements; j++)
+ if (strcmp(arg[i],elements[j]) == 0) break;
+ map[i-2] = j;
+ if (j == nelements) {
+ n = strlen(arg[i]) + 1;
+ elements[j] = new char[n];
+ strcpy(elements[j],arg[i]);
+ nelements++;
+ }
+ }
+
+ // read potential file and initialize potential parameters
+
+ read_file(arg[2]);
+
+ n = atom->ntypes;
+
+ // generate Coulomb 1/r energy look-up table
+
+ if (comm->me == 0 && screen) fprintf(screen,"Pair SMTBQ:\n");
+ if (comm->me == 0 && screen)
+ fprintf(screen," generating Coulomb integral lookup table ...\n");
+
+ tabqeq();
+
+ // ------------
+
+
+ if (comm->me == 0 && screen)
+ fprintf(screen," generating Second Moment integral lookup table ...\n");
+
+ tabsm();
+
+ // ------------
+
+ // clear setflag since coeff() called once with I,J = * *
+
+ 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");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+ ------------------------------------------------------------------------- */
+
+void PairSMTBQ::init_style()
+{
+ if (atom->tag_enable == 0)
+ error->all(FLERR,"Pair style SMTBQ requires atom IDs");
+ if (force->newton_pair == 0)
+ error->all(FLERR,"Pair style SMTBQ requires newton pair on");
+ if (!atom->q_flag)
+ error->all(FLERR,"Pair style SMTBQ requires atom attribute q");
+
+
+ // need a full neighbor list
+
+ int irequest = neighbor->request(this);
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->full = 1;
+
+ pgsize = neighbor->pgsize;
+ oneatom = neighbor->oneatom;
+ // if (maxpage == 0) add_pages();
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairSMTBQ::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+ return cutmax;
+}
+
+/* ----------------------------------------------------------------------
+ ---------------------------------------------------------------------- */
+
+void PairSMTBQ::read_file(char *file)
+{
+ int c, num_atom_types,i,k,m,test,j,verbose;
+ char **words;
+
+ memory->sfree(params);
+ params = NULL;
+ memory->sfree(intparams);
+ intparams = NULL;
+ nparams = 0;
+ maxparam = 0;
+ maxintparam = 0;
+
+ verbose = 1;
+ verbose = 0;
+
+ // open file on proc 0
+ FILE *fp;
+ fp = fopen( file, "r" );
+ if ( fp == NULL ) {
+ fprintf( stderr, "error opening the force filed file! terminating...\n" );
+ }
+
+
+ // read each line out of file, skipping blank lines or leading '#'
+ // store line of params if all 3 element tags are in element list
+
+ char *ptr;
+
+ ptr = (char*) malloc(sizeof(char)*MAXLINE);
+ words = (char**) malloc(sizeof(char*)*MAXTOKENS);
+ for (i=0; i < MAXTOKENS; i++)
+ words[i] = (char*) malloc(sizeof(char)*MAXTOKENS);
+
+
+ /* strip comment, skip line if blank */
+
+ if (verbose) printf ("\n");
+ fgets(ptr,MAXLINE,fp);
+ while (strchr(ptr,'#')) {
+ if (verbose) printf ("%s",ptr);
+ fgets(ptr,MAXLINE,fp);
+ }
+
+
+ // Nombre d'atome different dans la structure
+ // ===============================================
+ c = Tokenize( ptr, &words );
+ num_atom_types = atoi(words[1]);
+ if (verbose) printf (" %s %d\n", words[0], num_atom_types);
+
+ memory->create(intype,num_atom_types,num_atom_types,"pair:intype");
+
+ m = 0;
+ for (i = 0; i < num_atom_types; i++) {
+ for (j = 0; j < num_atom_types; j++) {
+ if (j < i) { intype[i][j] = intype[j][i];}
+ else { intype[i][j] = 0;
+ m = m + 1; }
+ if (verbose) printf ("i %d, j %d, intype %d - nb pair %d\n",i,j,intype[i][j],m);
+ }
+ }
+
+ // load up parameter settings and error check their values
+
+ if (nparams == maxparam) {
+ maxparam += DELTA;
+ params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
+ "pair:params");
+ maxintparam += m;
+ intparams = (Intparam *) memory->srealloc(intparams,(maxintparam+1)*sizeof(Intparam),
+ "pair:intparams");
+ }
+
+ for (i=0; i < num_atom_types; i++)
+ params[i].nom = (char*) malloc(sizeof(char)*3);
+
+ for (i=1; i <= maxintparam; i++)
+ intparams[i].typepot = (char*) malloc(sizeof(char)*15);
+
+ for (i=1; i <= maxintparam; i++)
+ intparams[i].mode = (char*) malloc(sizeof(char)*6);
+
+ QEqMode = (char*) malloc(sizeof(char)*18);
+ Bavard = (char*) malloc(sizeof(char)*5);
+ QInitMode = (char*) malloc(sizeof(char)*18);
+ writepot = (char*) malloc(sizeof(char)*5);
+ writeenerg = (char*) malloc(sizeof(char)*5);
+
+
+ // Little loop for ion's parameters
+ // ================================================
+ for (i=0; i %d = %s\n",words[1],i,params[i].nom);
+
+ for (j = 0; j %d = %s\n",words[2],j,params[j].nom);
+
+
+ if ( test == 1 ) {
+ if (verbose) printf ("========== fin des interaction ==========\n");
+ break ; }
+
+
+ intype[i][j] = m;
+ intype[j][i] = intype[i][j];
+ strcpy( intparams[m].typepot , words[3] );
+ intparams[m].intsm = 0;
+ if (verbose) printf (" itype %d jtype %d - intype %d\n",i,j,intype[i][j]);
+
+ if (strcmp(intparams[m].typepot,"second_moment") !=0 &&
+ strcmp(intparams[m].typepot,"buck") != 0 &&
+ strcmp(intparams[m].typepot,"buckPlusAttr") != 0) {
+ error->all(FLERR,"the potential other than second_moment or buckingham have not treated here\n");}
+
+
+ // On detemrine le type d'interaction
+ // -----------------------------------
+ if (strcmp(intparams[m].typepot,"second_moment") == 0) {
+ maxintsm += 1;
+ strcpy( intparams[m].mode , words[4] );
+ intparams[m].intsm = maxintsm;
+
+ if (strcmp(intparams[m].mode,"oxide") != 0 &&
+ strcmp(intparams[m].mode,"metal") != 0){
+ error->all(FLERR,"needs mode to second moment interaction : oxide or metal"); }
+
+// if (strcmp(intparams[m].mode,"oxide") == 0)
+// intparams[m].ncov = min((params[i].sto)*(params[i].n0),(params[j].sto)*(params[j].n0));
+
+ if (verbose) printf(" %s %s %s %s %s \n",words[0],words[1],words[2],
+ intparams[m].typepot,intparams[m].mode);
+
+ fgets( ptr, MAXLINE, fp);
+ c= Tokenize( ptr, &words );
+
+ intparams[m].a = atof(words[1]) ;
+ intparams[m].p = atof(words[2]) ;
+ intparams[m].ksi = atof(words[3]) ;
+ intparams[m].q = atof(words[4]) ;
+ if (verbose) printf (" %s %f %f %f %f\n",words[0],
+ intparams[m].a,intparams[m].p,intparams[m].ksi,intparams[m].q);
+
+ // Ligne 6 - rayon de coupure potentiel SM
+
+ fgets( ptr, MAXLINE, fp);
+ c= Tokenize( ptr, &words );
+
+ intparams[m].dc1 = atof(words[1]) ;
+ intparams[m].dc2 = atof(words[2]) ;
+ intparams[m].r0 = atof(words[3]) ;
+
+
+ if (strcmp(intparams[m].mode,"metal") == 0) {
+ if (verbose) printf (" %s %f %f %f\n",words[0],
+ intparams[m].dc1,intparams[m].dc2,intparams[m].r0);
+ } else {
+ if (verbose) printf (" %s %f %f %f\n",words[0],
+ intparams[m].dc1,intparams[m].dc2,intparams[m].r0);
+ }
+
+
+ } else if (strcmp(intparams[m].typepot,"buck") == 0) {
+
+ if (verbose) printf(" %s %s %s %s\n",words[0],words[1],words[2],
+ intparams[m].typepot);
+
+ fgets( ptr, MAXLINE, fp);
+ c= Tokenize( ptr, &words );
+
+ intparams[m].abuck = atof(words[1]) ; intparams[m].rhobuck = atof(words[2]) ;
+ if (verbose) printf (" %s %f %f\n",words[0],intparams[m].abuck,intparams[m].rhobuck);
+
+ }
+
+ else if (strcmp(intparams[m].typepot,"buckPlusAttr") == 0) {
+
+ if (verbose) printf(" %s %s %s %s\n",words[0],words[1],words[2],
+ intparams[m].typepot);
+
+ fgets( ptr, MAXLINE, fp);
+ c= Tokenize( ptr, &words );
+
+ intparams[m].abuck = atof(words[1]) ; intparams[m].rhobuck = atof(words[2]) ;
+ if (verbose) printf (" %s %f %f\n",words[0],intparams[m].abuck,intparams[m].rhobuck);
+
+
+ fgets( ptr, MAXLINE, fp);
+ c= Tokenize( ptr, &words );
+
+ intparams[m].aOO = atof(words[1]) ; intparams[m].bOO = atof(words[2]) ;
+ intparams[m].r1OO = atof(words[3]) ;intparams[m].r2OO = atof(words[4]) ;
+ if (verbose) printf (" %s %f %f %f %f \n",words[0],intparams[m].aOO,
+ intparams[m].bOO,intparams[m].r1OO,intparams[m].r2OO);
+
+ }
+ if (verbose) printf (" intsm %d \n",intparams[m].intsm);
+
+ } // for maxintparam
+
+
+ /* ====================================================================
+ tables Parameters
+ ==================================================================== */
+
+ // Ligne 9 - rayon de coupure Electrostatique
+ if (test == 0) {
+ fgets(ptr,MAXLINE,fp);
+ if (verbose) printf ("%s\n",ptr);
+
+ fgets( ptr, MAXLINE, fp);
+ }
+ c= Tokenize( ptr, &words );
+
+ for (i=0 ; iall(FLERR,"The QEq Mode is not known. QEq mode should be :\n Possible QEq modes | parameters\n QEqAll | no parameters\n QEqAllParallel | no parameters\n Surface | zlim (QEq only for z>zlim)\n BulkFromSlab | zlim1 zlim2 (QEq only for zlim1 %f Ang\n",
+ intparams[m].typepot,intparams[m].neig_cut);
+ }
+ }
+
+ //A adapter au STO
+ ncov = min((params[0].sto)*(params[0].n0),(params[1].sto)*(params[1].n0));
+
+ if (verbose) printf (" Parametre ncov = %f\n",ncov);
+ if (verbose) printf (" ********************************************* \n");
+
+
+
+}
+
+/* ----------------------------------------------------------------------
+ * COMPUTE
+ ---------------------------------------------------------------------- */
+
+void PairSMTBQ::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,m,gp;
+ int itag,jtag,itype,jtype;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+ double rsq,iq,jq,Eself,natom;
+ double ecovtot,ErepOO,ErepMO,Eion,Ecoh;
+ double **tmp,**tmpAll,*nmol;
+ double dq,dqcov;
+
+
+ int bavard;
+
+
+ if (atom->nmax > nmax) {
+ memory->destroy(ecov);
+ memory->destroy(potmad);
+ memory->destroy(potself);
+ memory->destroy(potcov);
+ memory->destroy(sbcov);
+ memory->destroy(coord);
+ memory->destroy(sbmet);
+ memory->destroy(chimet);
+ memory->destroy(flag_QEq);
+ memory->destroy(qf);
+ memory->destroy(q1);
+ memory->destroy(q2);
+ memory->destroy(tab_comm);
+
+ nmax = atom->nmax;
+
+ memory->create(ecov,nmax,"pair:ecov");
+ memory->create(potmad,nmax,"pair:potmad");
+ memory->create(potself,nmax,"pair:potself");
+ memory->create(potcov,nmax,"pair:potcov");
+ memory->create(sbcov,nmax,"pair:sbcov");
+ memory->create(coord,nmax,"pair:coord");
+ memory->create(sbmet,nmax,"pair:sbmet");
+ memory->create(chimet,nmax,"pair:chimet");
+ memory->create(flag_QEq,nmax,"pair:flag_QEq");
+ memory->create(qf,nmax,"pair:qf");
+ memory->create(q1,nmax,"pair:q1");
+ memory->create(q2,nmax,"pair:q2");
+ memory->create(tab_comm,nmax,"pair:tab_comm");
+ }
+
+
+ evdwl = ecoul = ecovtot = ErepOO = ErepMO = Eion = 0.0;
+ Eself = 0.0;
+
+ if (eflag || vflag) { ev_setup(eflag,vflag); }
+ else { evflag = vflag_fdotr = vflag_atom = 0; }
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *tag = atom->tag;
+ int *type = atom->type;
+ int newton_pair = force->newton_pair;
+ int nlocal = atom->nlocal;
+ int step = update->ntimestep;
+ struct timeval start, end;
+
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ if (step == 0 || (Qstep !=0 && fmod(double(step), double(Qstep)) == 0.0 )) Charge();
+// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+// this is necessary to get sbcov or sbmet table in order to caclulate the covalent or metal bonding
+ if (Qstep == 0 || fmod(double(step), double(Qstep)) != 0.0 ) QForce_charge(0);
+
+
+ // Charges Communication
+ // ----------------------
+ forward(q) ; reverse(q);
+
+ memory->create(nmol,nteam+1,"pair:nmol");
+ memory->create(tmp,nteam+1,7,"pair:tmp");
+ memory->create(tmpAll,nteam+1,7,"pair:tmpAll");
+
+
+ for (i=0; i ionic energy
+ 1 -> coulombian energy
+ 2 -> Electrosatic energy (ionic + Coulombian)
+ 3 -> Short int. Ox-Ox
+ 4 -> Short int. SMTB (repulsion)
+ 5 -> Covalent energy SMTB
+ 6 -> Somme des Q(i)²
+ ------------------------------------------------------------------------- */
+
+ // -----------
+ bavard = 0;
+ // -----------
+
+ /* -------------- N-body forces Calcul --------------- */
+
+ for (ii = 0; ii < inum; ii++) {
+// ===============================
+ i = ilist[ii];
+ itag = tag[i];
+ itype = map[type[i]];
+ iq = q[i];
+ gp = flag_QEq[i];
+
+ if (gp == 0 && itype > 0) natom += 1.0;
+ // if (gp == 0 && itype > 0) nmol[gp] += 1.0;
+
+
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+
+ // --- For atom i
+ tmp[gp][6] += iq*iq;
+
+
+// self energy, only on i atom
+// ---------------------------
+ Eself = self(¶ms[itype],iq);
+ tmp[gp][0] += Eself;
+ tmp[gp][2] += Eself;
+
+ if (evflag) ev_tally_full (i,0.0,2.0*Eself,0.0,0.0,0.0,0.0);
+
+
+// N-body energy of i
+// ---------------------
+ dq = fabs(params[itype].qform) - fabs(iq);
+ dqcov = dq*(2.0*ncov/params[itype].sto - dq);
+
+ ecov[i] = - sqrt(sbcov[i]*dqcov + sbmet[i]);
+ ecovtot += ecov[i];
+ tmp[gp][5] += ecov[i];
+
+ if (evflag) ev_tally_full(i,0.0,2.0*ecov[i],0.0,0.0,0.0,0.0);
+
+
+// Coulombian Interaction
+// -----------------------
+ evdwl = 2.0*Vself*iq*iq ;
+ tmp[gp][1] += Vself*iq*iq;
+ tmp[gp][2] += Vself*iq*iq;
+
+ if (evflag) ev_tally_full (i,0.0,evdwl,0.0,0.0,0.0,0.0);
+ evdwl = 0.0 ;
+
+
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+// ===============================
+ j = jlist[jj];
+ jtype = map[type[j]];
+ jtag = tag[j]; jq = q[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] < x[i][2]) 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;
+ }
+// .......................................................................
+
+
+ // # of interaction
+ // ----------------
+ m = intype[itype][jtype];
+
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+
+// ---------------------------------
+ if (sqrt(rsq) > cutmax) continue;
+// ---------------------------------
+
+
+ // Coulombian Energy
+ // ------------------
+ evdwl = 0.0 ; fpair = 0.0;
+ potqeq(i,j,iq,jq,rsq,fpair,eflag,evdwl);
+
+ tmp[gp][1] += evdwl;
+ tmp[gp][2] += evdwl;
+
+
+ // Coulombian Force
+ // -----------------
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+
+
+ if (evflag)
+ ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
+ evdwl = 0.0; fpair = 0.0 ;
+
+
+
+// ---------------------
+ if (m == 0) continue;
+// ---------------------
+
+// ----------------------------------------------
+ if ( strcmp(intparams[m].typepot,"buck") == 0 ||
+ strcmp(intparams[m].typepot,"buckPlusAttr") ==0 ) {
+// ----------------------------------------------
+
+ evdwl = 0.0; fpair =0.0;
+ rep_OO (&intparams[m],rsq,fpair,eflag,evdwl);
+ ErepOO += evdwl ;
+ tmp[gp][3] += evdwl;
+
+
+ // repulsion is pure two-body, sums up pair repulsive forces
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+
+
+ if (evflag)
+ ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
+ evdwl = 0.0; fpair = 0.0 ;
+
+ } // ----------------------------------- Rep O-O
+
+ if (strcmp(intparams[m].typepot,"buckPlusAttr") == 0 ) {
+// ----------------------------------------------
+
+ evdwl = 0.0; fpair =0.0;
+ Attr_OO (&intparams[m],rsq,fpair,eflag,evdwl);
+ ErepOO += evdwl ;
+ tmp[gp][3] += evdwl;
+
+
+ // repulsion is pure two-body, sums up pair repulsive forces
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+
+
+ if (evflag)
+ ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
+ evdwl = 0.0; fpair = 0.0 ;
+
+ } // ----------------------------------- Attr O-O
+
+
+// -----------------------------------------------------------------
+ if (strcmp(intparams[m].typepot,"second_moment") != 0 ) continue;
+// -----------------------------------------------------------------
+
+
+ if (sqrt(rsq) > intparams[m].dc2) continue;
+// -------------------------------------------
+
+// Repulsion : Energy + force
+// ----------------------------
+ evdwl = 0.0; fpair = 0.0 ;
+ repulsive(&intparams[m],rsq,i,j,fpair,eflag,evdwl);
+ ErepMO += evdwl;
+ tmp[gp][4] += 2.0*evdwl;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+
+
+ if (evflag)
+ ev_tally(i,j,nlocal,newton_pair,2.0*evdwl,0.0,fpair,delx,dely,delz);
+
+ evdwl = 0.0 ; fpair = 0.0;
+// ----- ----- ----- ----- ----- -----
+
+// Attraction : force
+// ------------------
+ fpair = 0.0;
+ f_att(&intparams[m], i, j, rsq, fpair) ;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+
+ if (evflag)
+ ev_tally(i,j,nlocal,newton_pair,0.0,0.0,fpair,delx,dely,delz);
+
+
+ } // --------------------------------- End j
+
+ } // ---------------------------------- End i
+
+
+ if (vflag_fdotr) virial_fdotr_compute();
+
+
+ for (i = 0; i < nteam+1; i++) {
+ MPI_Allreduce(tmp[i],tmpAll[i],7,MPI_DOUBLE,MPI_SUM,world);
+ }
+
+// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+ if (me == 0 && fmod(double(step), Nevery) == 0.0 && strcmp(writeenerg,"true") == 0) {
+
+ ofstream fichierE;
+
+ if (step == 0) { fichierE.open ("Energy_component.txt", ios::out | ios::trunc) ;}
+ else { fichierE.open ("Energy_component.txt", ios::out | ios::app) ;}
+
+ if (fichierE) fichierE<< setprecision(9) <pair->tail_flag,vflag_fdotr);
+ printf ("neighbor->includegroup %d\n",neighbor->includegroup);
+
+
+
+ for (gp=0; gpdestroy(nmol);
+ memory->destroy(tmp);
+ memory->destroy(tmpAll);
+
+}
+
+/* ----------------------------------------------------------------------
+ Partie Electrostatique
+ ----------------------------------------------------------------------*/
+
+double PairSMTBQ::self(Param *param, double qi)
+{
+ double self_tmp;
+ double s1=param->chi, s2=param->dj;
+
+ self_tmp = qi*(s1+0.5*qi*s2);
+
+ return self_tmp;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairSMTBQ::qfo_self(Param *param, double qi)
+{
+ double self_d;
+ double s1 = param->chi;
+ double s2 = param->dj;
+
+ self_d = 0.0 ;
+ self_d = s1+qi*s2;
+
+ return self_d;
+}
+
+/* ---------------------------------------------------------------------- */
+/* ---------------------------------------------------------------------- */
+
+void PairSMTBQ::tabqeq()
+{
+ int i,j,k,m,verbose;
+ int nntype;
+ double rc,s,r;
+ double alf;
+
+ int ii;
+ double za,zb,ra,rb,gam,dgam,dza,dzb,
+ d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2,na,nb;
+ double aCoeff,bCoeff,rcoupe,nang;
+
+ int n = atom->ntypes;
+ int nlocal = atom->nlocal;
+ int nghost = atom->nghost;
+ nmax = atom->nmax;
+
+ verbose = 1;
+ verbose = 0;
+
+ nntype = int((n+1)*n/2);
+
+ rc = cutmax ;
+ alf = 0.3 ;
+ // alf = 0.2 ;
+
+
+ if (verbose) printf ("kmax %d, ds %f, nmax %d\n",kmax,ds,nmax);
+ if (verbose) printf ("nlocal = %d, nghost = %d\n",nlocal,nghost);
+ if (verbose) printf ("nntypes %d, kmax %d, rc %f, n %d\n",nntype,kmax,rc,n);
+
+ // allocate arrays
+
+ memory->create(coultype,n,n,"pair:intype");
+ memory->create(potqn,kmax+5,"pair:potqn");
+ memory->create(dpotqn,kmax+5,"pair:dpotqn");
+ memory->create(fafb,kmax+5,nntype,"pair:fafb");
+ memory->create(dfafb,kmax+5,nntype,"pair:dfafb");
+ memory->create(fafbOxOxSurf,kmax+5,"pair:fafbOxOxSurf");
+ memory->create(dfafbOxOxSurf,kmax+5,"pair:dfafbOxOxSurf");
+ memory->create(fafbTiOxSurf,kmax+5,"pair:fafbTiOxSurf");
+ memory->create(dfafbTiOxSurf,kmax+5,"pair:dfafbTiOxSurf");
+
+ memory->create(fafbOxOxBB,kmax+5,"pair:fafbOxOxBB");
+ memory->create(dfafbOxOxBB,kmax+5,"pair:dfafbOxOxBB");
+ memory->create(fafbTiOxBB,kmax+5,"pair:fafbTiOxB");
+ memory->create(dfafbTiOxBB,kmax+5,"pair:dfafbTiOxBB");
+
+
+ memory->create(ecov,nmax,"pair:ecov");
+ memory->create(potmad,nmax,"pair:potmad");
+ memory->create(potself,nmax,"pair:potself");
+ memory->create(potcov,nmax,"pair:potcov");
+ memory->create(sbcov,nmax,"pair:sbcov");
+ memory->create(coord,nmax,"pair:coord");
+ memory->create(sbmet,nmax,"pair:sbmet");
+ memory->create(chimet,nmax,"pair:chimet");
+
+ // memory->create(nvsm,nmax,"pair:nvsm");
+ // memory->create(vsm,nmax,nmax,"pair:vsm");
+ memory->create(flag_QEq,nmax,"pair:flag_QEq");
+
+ memory->create(qf,nmax,"pair:qf");
+ memory->create(q1,nmax,"pair:q1");
+ memory->create(q2,nmax,"pair:q2");
+ memory->create(tab_comm,nmax,"pair:tab_comm");
+
+ memory->create(fct,31,"pair:fct");
+
+ // set interaction number: 0-0=0, 1-1=1, 0-1=1-0=2
+
+ m = 0; k = n;
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < n; j++) {
+ if (j == i) {
+ coultype[i][j] = m;
+ m += 1;
+ } else if (j != i && j > i) {
+ coultype[i][j] = k;
+ k += 1;
+ } else if (j != i && j < i) {
+ coultype[i][j] = coultype[j][i];
+ }
+ if (verbose) printf ("i %d, j %d, coultype %d\n",i,j,coultype[i][j]);
+ }
+ }
+
+ // -------- Tabqn --------
+
+ // -------------------
+ // Ouverture du fichier
+ // ofstream fichier("tabqeq.txt", ios::out | ios::trunc) ;
+ // -------------------
+
+ double pi,mu;
+ pi = 4.0*atan(1.0);
+
+ mu = erfc(alf*rc)/rc ;
+
+ //if (fichier) fichier <<" r - potqn " <rcoupe) nang = rcoupe - rc ;
+ bCoeff = (2*dij+ddij*nang)/(dij*nang);
+ aCoeff = dij*exp(-bCoeff*rc) /pow(nang,2);
+ }
+ if (r > rc) {dij = aCoeff *pow((r- rc-nang),2) *exp(bCoeff*r);
+ ddij = aCoeff*(r- rc-nang) *(2+bCoeff*(r-rc-nang))*exp(bCoeff*r);
+ }
+
+
+ if (r > (rc+nang)) {dij = 0.0 ; ddij = 0.0;}
+
+ fafb[k][m] = potqn[k] - dij ;
+ if (k == 1) fafb[0][m] = fafb[k][m] ;
+
+ dfafb[k][m] = dpotqn[k] - ddij/r ;
+ }
+
+ // Make the table fafbOxOxSurf
+ rc = cutmax;
+ if(strcmp(params[i].nom,"O")==0 || strcmp(params[j].nom,"O")==0){
+ if(strcmp(params[i].nom,"O")==0) {
+ ra = ROxSurf;
+ za = (2.0*params[i].ne + 1.0)/(4.0*ra);}
+ if(strcmp(params[j].nom,"O")==0) {
+ rb = ROxSurf;
+ zb = (2.0*params[j].ne + 1.0)/(4.0*rb); }
+
+
+ ii = 0 ; nang =cang= 5.0 ;
+ // --------------------------
+ for (k = 0; k < kmax+5; k++)
+ // --------------------------
+ {
+ gam = dgam = dza = dzb = d2zaa = d2zab =
+ d2zbb = d2zra = d2zrb = d2gamr2 = 0.0 ;
+ dij = 0.0 ;
+
+ s = double(k)*ds ; r = sqrt(s) ;
+ if (k==0) r=10e-30;
+
+ gammas(na,nb,za,zb,r,gam,dgam,dza,dzb,
+ d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
+
+ // --- Jij
+
+ dij = 14.4 * (1.0/r - double(gam));
+ ddij = 14.4 * (-1.0/(r*r) - double(dgam)) ;
+
+ if (dij < 0.01 && ii==0)
+ {
+ ii=2;
+ if (ii==2) if (verbose) printf ("rc : %f\n",r);
+ rc = r ; ii=1 ;
+ if ((rc+nang)>rcoupe) nang = rcoupe - rc ;
+ bCoeff = (2*dij+ddij*nang)/(dij*nang);
+ aCoeff = dij*exp(-bCoeff*rc) /pow(nang,2);
+ }
+ if (r > rc) {dij = aCoeff *pow((r- rc-nang),2) *exp(bCoeff*r);
+ ddij = aCoeff*(r- rc-nang) *(2+bCoeff*(r-rc-nang))*exp(bCoeff*r);
+ }
+
+
+ if (r > (rc+nang)) {dij = 0.0 ; ddij = 0.0;}
+
+ if(strcmp(params[i].nom,"O")==0 && strcmp(params[j].nom,"O")==0){
+ fafbOxOxSurf[k] = potqn[k] - dij ;
+ if (k == 1) fafbOxOxSurf[0] = fafbOxOxSurf[k] ;
+
+ dfafbOxOxSurf[k] = dpotqn[k] - ddij/r ;
+ }
+ else { fafbTiOxSurf[k] = potqn[k] - dij ;
+ if (k == 1) fafbTiOxSurf[0] = fafbTiOxSurf[k] ;
+
+ dfafbTiOxSurf[k] = dpotqn[k] - ddij/r ;}
+
+ }
+
+
+ } //for k
+ //end of make the table fafbOxOxSurf
+
+ // Makes the table fafbOxOxBB
+ rc = cutmax;
+ if(strcmp(params[i].nom,"O")==0 || strcmp(params[j].nom,"O")==0){
+ if(strcmp(params[i].nom,"O")==0) {
+ ra = ROxBB;
+ za = (2.0*params[i].ne + 1.0)/(4.0*ra);}
+ if(strcmp(params[j].nom,"O")==0) {
+ rb = ROxBB;
+ zb = (2.0*params[j].ne + 1.0)/(4.0*rb); }
+
+
+ ii = 0 ; nang =cang= 5.0 ;
+ // --------------------------
+ for (k = 0; k < kmax+5; k++)
+ // --------------------------
+ {
+ gam = dgam = dza = dzb = d2zaa = d2zab =
+ d2zbb = d2zra = d2zrb = d2gamr2 = 0.0 ;
+ dij = 0.0 ;
+
+ s = double(k)*ds ; r = sqrt(s) ;
+ if (k==0) r=10e-30;
+
+ gammas(na,nb,za,zb,r,gam,dgam,dza,dzb,
+ d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
+
+ // --- Jij
+
+ dij = 14.4 * (1.0/r - double(gam));
+ ddij = 14.4 * (-1.0/(r*r) - double(dgam)) ;
+
+ if (dij < 0.01 && ii==0) {
+ ii=2;
+ if (ii==2) if (verbose) printf ("rc : %f\n",r);
+ rc = r ; ii=1 ;
+ if ((rc+nang)>rcoupe) nang = rcoupe - rc ;
+ bCoeff = (2*dij+ddij*nang)/(dij*nang);
+ aCoeff = dij*exp(-bCoeff*rc) /pow(nang,2);
+ }
+ if (r > rc) {
+ dij = aCoeff *pow((r- rc-nang),2) *exp(bCoeff*r);
+ ddij = aCoeff*(r- rc-nang) *(2+bCoeff*(r-rc-nang))*exp(bCoeff*r);
+ }
+
+
+ if (r > (rc+nang)) {dij = 0.0 ; ddij = 0.0;}
+
+ if(strcmp(params[i].nom,"O")==0 && strcmp(params[j].nom,"O")==0){
+ fafbOxOxBB[k] = potqn[k] - dij ;
+ if (k == 1) fafbOxOxBB[0] = fafbOxOxBB[k] ;
+ dfafbOxOxBB[k] = dpotqn[k] - ddij/r ; }
+ else { fafbTiOxBB[k] = potqn[k] - dij ;
+ if (k == 1) fafbTiOxBB[0] = fafbTiOxBB[k] ;
+ dfafbTiOxBB[k] = dpotqn[k] - ddij/r ;
+ }
+ }
+
+
+
+ } //for k
+ //end of make the table fafbOxOxBB
+
+
+
+ }
+ } //for i,j
+
+ //if (fichier) fichier.close() ;
+
+}
+
+/* ---------------------------------------------------------------------*/
+
+void PairSMTBQ::potqeq(int i, int j, double qi, double qj, double rsq,
+ double &fforce, int eflag, double &eng)
+{
+
+ /* ===================================================================
+ Coulombian energy calcul between i and j atoms
+ with fafb table make in sm_table().
+ fafb[i][j] : i is the table's step (r)
+ j is the interaction's # (in intype[itype][jtype])
+ dfafb is derivate energy (force)
+ ==================================================================== */
+
+ int itype,jtype,l,m;
+ double r,t1,t2,sds,xi,engBulk,engSurf,fforceBulk,fforceSurf,dcoordloc,dcoupureloc;
+ double engBB,fforceBB, dIntfcoup2loc,iCoord,jCoord,iIntfCoup2,jIntfCoup2;
+
+ int *type = atom->type;
+ // int n = atom->ntypes;
+
+ itype = map[type[i]];
+ jtype = map[type[j]];
+ m = coultype[itype][jtype];
+
+ r = rsq;
+ sds = r/ds ; l = int(sds) ;
+ xi = sds - double(l) ;
+
+
+ iCoord=coord[i];
+ jCoord=coord[j];
+ iIntfCoup2= Intfcoup2(iCoord,coordOxBulk,0.15);
+ jIntfCoup2= Intfcoup2(jCoord,coordOxBulk,0.15);
+
+ // ---- Energies Interpolation
+
+ t1 = fafb[l][m] + (fafb[l+1][m] - fafb[l][m])*xi;
+ t2 = fafb[l+1][m] + (fafb[l+2][m] - fafb[l+1][m])*(xi-1.0);
+
+ engBulk = qi*qj*(t1 + (t2 - t1)*xi/2.0);
+ eng=engBulk;
+
+
+ // ---- Forces Interpolation
+
+ t1 = dfafb[l][m] + (dfafb[l+1][m] - dfafb[l][m])*xi;
+ t2 = dfafb[l+1][m] + (dfafb[l+2][m] - dfafb[l+1][m])*(xi-1);
+
+
+ fforce = - qi*qj*(t1 + (t2 - t1)*xi/2.0) ;
+
+
+ if(strcmp(params[itype].nom,"O")==0 || strcmp(params[jtype].nom,"O")==0){
+
+ if(strcmp(params[itype].nom,"O")==0 && strcmp(params[jtype].nom,"O")==0){
+ // between two oxygens
+
+ t1 = fafbOxOxSurf[l] + (fafbOxOxSurf[l+1] - fafbOxOxSurf[l])*xi;
+ t2 = fafbOxOxSurf[l+1] + (fafbOxOxSurf[l+2] - fafbOxOxSurf[l+1])*(xi-1.0);
+ engSurf = qi*qj*(t1 + (t2 - t1)*xi/2.0);
+
+ t1 = fafbOxOxBB[l] + (fafbOxOxBB[l+1] - fafbOxOxBB[l])*xi;
+ t2 = fafbOxOxBB[l+1] + (fafbOxOxBB[l+2] - fafbOxOxBB[l+1])*(xi-1.0);
+ engBB = qi*qj*(t1 + (t2 - t1)*xi/2.0);
+
+ eng= engBulk + (iCoord+jCoord-2*coordOxBulk)/(2*(coordOxBB-coordOxBulk)) *(engBB-engBulk)
+ + (iIntfCoup2+jIntfCoup2)*((engBulk-engSurf)/(2*(coordOxBulk-coordOxSurf))
+ - (engBB-engBulk)/(2*(coordOxBB-coordOxBulk))) ;
+
+
+ // ---- Interpolation des forces
+
+ fforceBulk=fforce;
+ t1 = dfafbOxOxSurf[l] + (dfafbOxOxSurf[l+1] - dfafbOxOxSurf[l])*xi;
+ t2 = dfafbOxOxSurf[l+1] + (dfafbOxOxSurf[l+2] - dfafbOxOxSurf[l+1])*(xi-1);
+ fforceSurf = - qi*qj*(t1 + (t2 - t1)*xi/2.0) ;
+
+ t1 = dfafbOxOxBB[l] + (dfafbOxOxBB[l+1] - dfafbOxOxBB[l])*xi;
+ t2 = dfafbOxOxBB[l+1] + (dfafbOxOxBB[l+2] - dfafbOxOxBB[l+1])*(xi-1);
+ fforceBB = - qi*qj*(t1 + (t2 - t1)*xi/2.0) ;
+
+ fforce= fforceBulk + (iCoord+jCoord-2*coordOxBulk)/(2*(coordOxBB-coordOxBulk))*(fforceBB-fforceBulk)
+ + (iIntfCoup2+jIntfCoup2)*((fforceBulk-fforceSurf)/(2*(coordOxBulk-coordOxSurf))
+ - (fforceBB-fforceBulk)/(2*(coordOxBB-coordOxBulk))) ;
+
+ }
+ else{ // between metal and oxygen
+
+ t1 = fafbTiOxSurf[l] + (fafbTiOxSurf[l+1] - fafbTiOxSurf[l])*xi;
+ t2 = fafbTiOxSurf[l+1] + (fafbTiOxSurf[l+2] - fafbTiOxSurf[l+1])*(xi-1.0);
+ engSurf = qi*qj*(t1 + (t2 - t1)*xi/2.0);
+ t1 = fafbTiOxBB[l] + (fafbTiOxBB[l+1] - fafbTiOxBB[l])*xi;
+ t2 = fafbTiOxBB[l+1] + (fafbTiOxBB[l+2] - fafbTiOxBB[l+1])*(xi-1.0);
+ engBB = qi*qj*(t1 + (t2 - t1)*xi/2.0);
+
+ if(strcmp(params[jtype].nom,"O")==0) //the atom j is an oxygen
+ { iIntfCoup2=jIntfCoup2;
+ iCoord=jCoord; }
+
+ eng = engBulk + (engBulk-engSurf)/(coordOxBulk-coordOxSurf) * iIntfCoup2
+ + (engBB-engBulk)/(coordOxBB-coordOxBulk) * (iCoord-coordOxBulk-iIntfCoup2);
+
+
+ // ---- Forces Interpolation
+
+ fforceBulk=fforce;
+ t1 = dfafbTiOxSurf[l] + (dfafbTiOxSurf[l+1] - dfafbTiOxSurf[l])*xi;
+ t2 = dfafbTiOxSurf[l+1] + (dfafbTiOxSurf[l+2] - dfafbTiOxSurf[l+1])*(xi-1);
+ fforceSurf = - qi*qj*(t1 + (t2 - t1)*xi/2.0) ;
+
+ t1 = dfafbTiOxBB[l] + (dfafbTiOxBB[l+1] - dfafbTiOxBB[l])*xi;
+ t2 = dfafbTiOxBB[l+1] + (dfafbTiOxBB[l+2] - dfafbTiOxBB[l+1])*(xi-1);
+ fforceBB = - qi*qj*(t1 + (t2 - t1)*xi/2.0) ;
+
+ dcoordloc = fcoupured(sqrt(r),r1Coord,r2Coord) ;
+
+
+ dcoupureloc = fcoupured(iCoord,coordOxSurf,coordOxBulk) ;
+ dIntfcoup2loc= fcoup2(iCoord,coordOxBulk,0.15)*dcoupureloc ;
+ fforce = fforceBulk + 1/(coordOxBulk-coordOxSurf) * ((fforceBulk-fforceSurf)* iIntfCoup2
+ - (engBulk-engSurf) *dIntfcoup2loc)
+ + 1/(coordOxBB-coordOxBulk) * ((fforceBB-fforceBulk)*(iCoord-coordOxBulk- iIntfCoup2)
+ - (engBB-engBulk) *(dcoordloc-dIntfcoup2loc));
+
+
+ }
+
+
+
+ }
+
+
+}
+
+/* -------------------------------------------------------------------- */
+
+void PairSMTBQ::pot_ES (int i, int j, double rsq, double &eng)
+{
+
+ /* ===================================================================
+ Coulombian potentiel energy calcul between i and j atoms
+ with fafb table make in sm_table().
+ fafb[i][j] : i is the table's step (r)
+ j is the interaction's # (in intype[itype][jtype])
+ dfafb is derivate energy (force)
+ ==================================================================== */
+
+ int itype,jtype,l,m;
+ double r,t1,t2,sds,xi,engBulk,engSurf,dcoordloc,dcoupureloc;
+ double engBB, dIntfcoup2loc,iCoord,jCoord,iIntfCoup2,jIntfCoup2;
+
+ int *type = atom->type;
+ // int n = atom->ntypes;
+
+ itype = map[type[i]];
+ jtype = map[type[j]];
+ m = coultype[itype][jtype];
+
+ r = rsq;
+ sds = r/ds ; l = int(sds) ;
+ xi = sds - double(l) ;
+
+
+ iCoord=coord[i];
+ jCoord=coord[j];
+ iIntfCoup2= Intfcoup2(iCoord,coordOxBulk,0.15);
+ jIntfCoup2= Intfcoup2(jCoord,coordOxBulk,0.15);
+
+ // ---- Energies Interpolation
+
+ t1 = fafb[l][m] + (fafb[l+1][m] - fafb[l][m])*xi;
+ t2 = fafb[l+1][m] + (fafb[l+2][m] - fafb[l+1][m])*(xi-1.0);
+
+
+ eng = (t1 + (t2 - t1)*xi/2.0);
+ engBulk=eng;
+
+
+ if(itype==0 || jtype==0){
+
+ if(itype==0 && jtype==0){ // between two oxygens
+
+ t1 = fafbOxOxSurf[l] + (fafbOxOxSurf[l+1] - fafbOxOxSurf[l])*xi;
+ t2 = fafbOxOxSurf[l+1] + (fafbOxOxSurf[l+2] - fafbOxOxSurf[l+1])*(xi-1.0);
+ engSurf = (t1 + (t2 - t1)*xi/2.0);
+
+ t1 = fafbOxOxBB[l] + (fafbOxOxBB[l+1] - fafbOxOxBB[l])*xi;
+ t2 = fafbOxOxBB[l+1] + (fafbOxOxBB[l+2] - fafbOxOxBB[l+1])*(xi-1.0);
+ engBB = (t1 + (t2 - t1)*xi/2.0);
+
+ eng= engBulk + (iCoord+jCoord-2*coordOxBulk)/(2*(coordOxBB-coordOxBulk))*(engBB-engBulk)
+ + (iIntfCoup2+jIntfCoup2)*((engBulk-engSurf)/(2*(coordOxBulk-coordOxSurf))
+ - (engBB-engBulk)/(2*(coordOxBB-coordOxBulk))) ;
+
+ } else { // between metal and oxygen
+
+ t1 = fafbTiOxSurf[l] + (fafbTiOxSurf[l+1] - fafbTiOxSurf[l])*xi;
+ t2 = fafbTiOxSurf[l+1] + (fafbTiOxSurf[l+2] - fafbTiOxSurf[l+1])*(xi-1.0);
+ engSurf = (t1 + (t2 - t1)*xi/2.0);
+
+ t1 = fafbTiOxBB[l] + (fafbTiOxBB[l+1] - fafbTiOxBB[l])*xi;
+ t2 = fafbTiOxBB[l+1] + (fafbTiOxBB[l+2] - fafbTiOxBB[l+1])*(xi-1.0);
+ engBB = (t1 + (t2 - t1)*xi/2.0);
+
+ if (jtype==0) { //the atom j is an oxygen
+ iIntfCoup2=jIntfCoup2;
+ iCoord=jCoord;
+ }
+
+ eng = engBulk + (engBulk-engSurf)/(coordOxBulk-coordOxSurf)*iIntfCoup2
+ + (engBB-engBulk)/(coordOxBB-coordOxBulk) * (iCoord-coordOxBulk-iIntfCoup2);
+
+
+ }
+
+
+
+ }
+
+
+}
+
+/* -------------------------------------------------------------------- */
+
+void PairSMTBQ::pot_ES2 (int i, int j, double rsq, double &pot)
+{
+ int l,m,itype,jtype;
+ double sds,xi,t1,t2,r;
+
+ int *type = atom->type;
+
+
+ if (sqrt(rsq) > cutmax) return ;
+
+ itype = map[type[i]];
+ jtype = map[type[j]];
+ m = coultype[itype][jtype];
+
+ r = rsq ;
+ sds = r/ds ; l = int(sds) ;
+ xi = sds - double(l) ;
+
+ // ---- Energies Interpolation
+
+ t1 = fafb[l][m] + (fafb[l+1][m] - fafb[l][m])*xi;
+ t2 = fafb[l+1][m] + (fafb[l+2][m] - fafb[l+1][m])*(xi-1.0);
+
+ pot = (t1 + (t2 - t1)*xi/2.0) ;
+
+}
+
+/* --------------------------------------------------------------------
+ Oxygen-Oxygen Interaction
+ -------------------------------------------------------------------- */
+
+void PairSMTBQ::rep_OO(Intparam *intparam, double rsq, double &fforce,
+ int eflag, double &eng)
+{
+ double r,tmp_exp,tmp;
+ double A = intparam->abuck ;
+ double rho = intparam->rhobuck ;
+
+ r = sqrt(rsq);
+ tmp = - r/rho ;
+ tmp_exp = exp( tmp );
+
+ eng = A * tmp_exp ;
+
+ fforce = A/rho * tmp_exp / r ; //( - )
+
+}
+
+
+void PairSMTBQ::Attr_OO(Intparam *intparam, double rsq, double &fforce,
+ int eflag, double &eng)
+{
+ double r,tmp_exp,tmp;
+ double aOO = intparam->aOO ;
+ double bOO = intparam->bOO ;
+ double r1OO = intparam->r1OO ;
+ double r2OO = intparam->r2OO ;
+
+ r = sqrt(rsq);
+ tmp_exp= exp( bOO* r);
+ eng = aOO * tmp_exp* fcoupure(r,r1OO,r2OO);
+
+ fforce = - (aOO*bOO * tmp_exp * fcoupure(r,r1OO,r2OO)+ aOO*tmp_exp *fcoupured(r,r1OO,r2OO))/ r ; //( - )
+
+}
+
+
+/* ----------------------------------------------------------------------
+ covalente Interaction
+ ----------------------------------------------------------------------*/
+
+
+void PairSMTBQ::tabsm()
+{
+ int k,m;
+ double s,r,tmpb,tmpr,fcv,fcdv;
+
+ memory->create(tabsmb,kmax,maxintsm+1,"pair:tabsmb");
+ memory->create(tabsmr,kmax,maxintsm+1,"pair:tabsmr");
+ memory->create(dtabsmb,kmax,maxintsm+1,"pair:dtabsmb");
+ memory->create(dtabsmr,kmax,maxintsm+1,"pair:dtabsmr");
+
+
+ for (m = 0; m <= maxintparam; m++) {
+
+ if (intparams[m].intsm == 0) continue;
+
+ double rc1 = intparams[m].dc1;
+ double rc2 = intparams[m].dc2;
+ double A = intparams[m].a;
+ double p = intparams[m].p;
+ double Ksi = intparams[m].ksi;
+ double q = intparams[m].q;
+ double rzero = intparams[m].r0;
+ int sm = intparams[m].intsm;
+
+
+ for (k=0; k < kmax; k++)
+ {
+ s = double(k)*ds ; r = sqrt(s);
+ if (k==0) r=10e-30;
+ tmpb = exp( -2.0*q*(r/rzero - 1.0));
+ tmpr = exp( -p*(r/rzero - 1.0));
+
+ if (r <= rc1)
+ {
+
+ // -- Energy
+ tabsmb[k][sm] = Ksi*Ksi * tmpb ;
+ tabsmr[k][sm] = A * tmpr ;
+
+ // -- Force
+ /* dtabsmb ne correspond pas vraiment a une force puisqu'il y a le /r
+ (on a donc une unite force/distance). Le programme multiplie ensuite
+ (dans le PairSMTBQ::compute ) dtabsmb par la projection du vecteur r
+ sur un axe x (ou y ou z) pour determiner la composante de la force selon
+ cette direction. Donc tout est ok au final. */
+
+ dtabsmb[k][sm] = - 2.0 *Ksi*Ksi* q/rzero * tmpb /r;
+ dtabsmr[k][sm] = - A * p/rzero * tmpr/r ;
+
+ } // if
+
+ else if (r > rc1 && r <= rc2)
+ {
+
+ // -- Energie
+ fcv = fcoupure(r,intparams[sm].dc1,intparams[sm].dc2);
+ tabsmb[k][sm] = fcv* Ksi*Ksi * tmpb ;
+ tabsmr[k][sm] = fcv* A * tmpr ;
+
+ // -- Force
+ /* dtabsmb ne correspond pas vraiment a une force puisqu'il y a le /r
+ (on a donc une unite force/distance). Le programme multiplie ensuite
+ (dans le PairSMTBQ::compute ) d tabsmb par la projection du vecteur
+ r sur un axe x (ou y ou z) pour determiner la composante de la force
+ selon cette direction. Donc tout est ok au final. */
+
+ fcdv = fcoupured(r,intparams[sm].dc1,intparams[sm].dc2);
+ dtabsmb[k][sm] = (fcv*( - 2.0 *Ksi*Ksi* q/rzero * tmpb )+fcdv* Ksi*Ksi * tmpb )/r ;
+ dtabsmr[k][sm] = (fcv*( - A * p/rzero * tmpr )+fcdv*A * tmpr )/r ;
+
+ }
+
+ else
+ {
+
+ // -- Energie
+ tabsmb[k][sm] = 0.0;
+ tabsmr[k][sm] = 0.0;
+
+ // -- Force
+ dtabsmb[k][sm] = 0.0;
+ dtabsmr[k][sm] = 0.0;
+
+ }
+
+
+
+ } // for kmax
+
+
+ } // for maxintparam
+
+}
+
+
+
+
+
+/* -------------------------------------------------------------- */
+
+void PairSMTBQ::repulsive(Intparam *intparam, double rsq, int i, int j,
+ double &fforce, int eflag, double &eng)
+{
+
+ /* ================================================
+ rsq : square of ij distance
+ fforce : repulsion force
+ eng : repulsion energy
+ eflag : Si oui ou non on calcule l'energie
+ =================================================*/
+
+ int l;
+ double r,sds,xi,t1,t2,dt1,dt2,sweet,iq,jq;
+
+ double rrcs = intparam->dc2;
+ int sm = intparam->intsm;
+
+ // printf ("On rentre dans repulsive\n");
+
+
+ r = rsq;
+ if (sqrt(r) > rrcs) return ;
+
+ sds = r/ds ; l = int(sds) ;
+ xi = sds - double(l) ;
+
+ t1 = tabsmr[l][sm] + (tabsmr[l+1][sm] - tabsmr[l][sm])*xi ;
+ t2 = tabsmr[l+1][sm] + (tabsmr[l+2][sm] - tabsmr[l+1][sm])*(xi-1.0) ;
+
+ dt1 = dtabsmr[l][sm] + (dtabsmr[l+1][sm] - dtabsmr[l][sm])*xi ;
+ dt2 = dtabsmr[l+1][sm] + (dtabsmr[l+2][sm] - dtabsmr[l+1][sm])*(xi-1.0) ;
+
+ if (strcmp(intparam->mode,"oxide") == 0)
+ {
+ fforce = - 2.0*(dt1 + (dt2 - dt1)*xi/2.0);
+ eng = (t1 + (t2 - t1)*xi/2.0) ;
+ }
+ else if (strcmp(intparam->mode,"metal") == 0)
+ {
+ sweet = 1.0;
+ fforce = - 2.0*(dt1 + (dt2 - dt1)*xi/2.0) * sweet ;
+ eng = (t1 + (t2 - t1)*xi/2.0) * sweet ;
+ }
+
+}
+
+
+/* --------------------------------------------------------------------------------- */
+
+
+void PairSMTBQ::attractive(Intparam *intparam, double rsq,
+ int eflag, int i, double iq, int j, double jq)
+{
+ int itype,l;
+ double r,t1,t2,xi,sds;
+ double dqcov,sweet,dq,mu;
+
+ double rrcs = intparam->dc2;
+ int *type = atom->type;
+ int sm = intparam->intsm;
+
+ itype = map[type[i]];
+
+ r = rsq;
+ if (sqrt(r) > rrcs) return ;
+
+
+ sds = r/ds ; l = int(sds) ;
+ xi = sds - double(l) ;
+
+ t1 = tabsmb[l][sm] + (tabsmb[l+1][sm] - tabsmb[l][sm])*xi ;
+ t2 = tabsmb[l+1][sm] + (tabsmb[l+2][sm] - tabsmb[l+1][sm])*(xi-1.0) ;
+
+
+
+ if (strcmp(intparam->mode,"oxide") == 0) {
+ mu = 0.5*(sqrt(params[1].sto) + sqrt(params[0].sto));
+
+// dq = fabs(params[itype].qform) - fabs(iq);
+// dqcov = dq*(2.0*ncov/params[itype].sto - dq);
+
+ sbcov[i] += (t1 + (t2 - t1)*xi/2.0) *params[itype].sto*mu*mu;
+
+// if (i < 10) printf ("i %d, iq %f sbcov %f \n",i,iq,sbcov[i]);
+
+ if (sqrt(r)mode,"metal") == 0) {
+ sweet = 1.0;
+ sbmet[i] += (t1 + (t2 - t1)*xi/2.0) * sweet ;
+ }
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSMTBQ::f_att(Intparam *intparam, int i, int j,double rsq, double &fforce)
+{
+ int itype,jtype,l;
+ int *type = atom->type;
+
+ double r,sds,xi,dt1,dt2,dq,dqcovi,dqcovj;
+ double fcov_ij,fcov_ji,sweet,iq,jq,mu;
+
+ int sm = intparam->intsm;
+ double *q = atom->q;
+
+ itype = map[type[i]];
+ jtype = map[type[j]];
+ iq = q[i] ; jq = q[j];
+
+ r = rsq;
+
+ sds = r/ds ; l = int(sds) ;
+ xi = sds - double(l) ;
+
+ dt1 = dtabsmb[l][sm] + (dtabsmb[l+1][sm] - dtabsmb[l][sm])*xi ;
+ dt2 = dtabsmb[l+1][sm] + (dtabsmb[l+2][sm] - dtabsmb[l+1][sm])*(xi-1.0) ;
+
+ dq = fabs(params[itype].qform) - fabs(iq);
+ dqcovi = dq*(2.0*ncov/params[itype].sto - dq);
+
+ dq = fabs(params[jtype].qform) - fabs(jq);
+ dqcovj = dq*(2.0*ncov/params[jtype].sto - dq);
+
+ if (strcmp(intparam->mode,"oxide") == 0) {
+//------------------------------------------
+ mu = 0.5*(sqrt(params[1].sto) + sqrt(params[0].sto));
+ fcov_ij = (dt1 + (dt2 - dt1)*xi/2.0) * dqcovi *params[itype].sto*mu*mu;
+ fcov_ji = (dt1 + (dt2 - dt1)*xi/2.0) * dqcovj *params[jtype].sto*mu*mu;
+
+ fforce = 0.5 * ( fcov_ij/sqrt(sbcov[i]*dqcovi + sbmet[i])
+ + fcov_ji/sqrt(sbcov[j]*dqcovj + sbmet[j]) ) ;
+ }
+
+ else if (strcmp(intparam->mode,"metal") == 0) {
+//-----------------------------------------------
+ sweet = 1.0;
+ fcov_ij = (dt1 + (dt2 - dt1)*xi/2.0) * sweet ;
+
+ fforce = 0.5 * fcov_ij*( 1.0/sqrt(sbcov[i]*dqcovi + sbmet[i])
+ + 1.0/sqrt(sbcov[j]*dqcovj + sbmet[j]) ) ;
+ }
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSMTBQ::pot_COV(Param *param, int i, double &qforce)
+{
+ double iq,dq,DQ,sign;
+
+ double *q = atom->q;
+ double qform = param->qform;
+ double sto = param->sto;
+
+ sign = qform / fabs(qform);
+ iq = q[i];
+
+ dq = fabs(qform) - fabs(iq);
+ DQ = dq*(2.0*ncov/sto - dq);
+
+ if (fabs(iq) < 1.0e-7 || fabs(sbcov[i]) < 1.0e-7) {
+ qforce = 0.0; }
+ else {
+ qforce = sign*sbcov[i]/sqrt(sbcov[i]*DQ + sbmet[i])*(ncov/sto - dq) ;
+ }
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairSMTBQ::potmet(Intparam *intparam, double rsq,
+ int i, double iq, int j, double jq)
+{
+ int l,itype,jtype;
+ int *type = atom->type;
+ double chi,sds,xi,t1,t2,r,dsweet,dq,dqcovi,dqcovj;
+
+ int sm = intparam->intsm;
+ itype = map[type[i]];
+ jtype = map[type[j]];
+
+ r = rsq;
+ sds = r/ds ; l = int(sds) ;
+ xi = sds - double(l) ;
+
+ t1 = tabsmb[l][sm] + (tabsmb[l+1][sm] - tabsmb[l][sm])*xi ;
+ t2 = tabsmb[l+1][sm] + (tabsmb[l+2][sm] - tabsmb[l+1][sm])*(xi-1.0) ;
+
+ dq = fabs(params[itype].qform) - fabs(iq);
+ dqcovi = dq*(2.0*ncov/params[itype].sto - dq);
+
+ dq = fabs(params[jtype].qform) - fabs(jq);
+ dqcovj = dq*(2.0*ncov/params[jtype].sto - dq);
+
+ dsweet = 0.0;
+ chi = (t1 + (t2 - t1)*xi/2.0) * dsweet *( 1.0/(2.0*sqrt(sbcov[i]*dqcovi+sbmet[i]))
+ + 1.0/(2.0*sqrt(sbcov[j]*dqcovj+sbmet[j])) );
+
+ return chi;
+}
+
+
+/* ----------------------------------------------------------------------
+ Cutting Function
+ ----------------------------------------------------------------------*/
+
+
+/* -------------------------------------------------------------------- */
+
+double PairSMTBQ::fcoupure(double r, double rep_dc1, double rep_dc2)
+{
+ double ddc,a3,a4,a5,x;
+
+ if (r rep_dc2)
+ {return 0;}
+ else
+ {
+
+ ddc = rep_dc2 - rep_dc1;
+ x = r - rep_dc2;
+
+ a3 = -10/(ddc*ddc*ddc);
+ a4 = -15/(ddc*ddc*ddc*ddc);
+ a5 = -6/(ddc*ddc*ddc*ddc*ddc);
+
+ return x*x*x*(a3 + x*(a4 + x*a5));}
+}
+
+/* ----------------------------------------------------------------------
+ Derivate of cutting function
+ ----------------------------------------------------------------------*/
+
+
+/* ----------------------------------------------------------------------- */
+
+double PairSMTBQ::fcoupured(double r, double rep_dc1, double rep_dc2)
+{
+
+ double ddc,a3,a4,a5,x;
+
+ if ( r>rep_dc1 && r x+delta)
+ {return 0;}
+ else
+ {
+ dc = c - x-delta;
+ return dc*dc*(3*delta+dc)/(4*delta*delta*delta);
+ }
+}
+
+/* ----------------------------------------------------------------------
+ Primitive of cutting function for derive (force)
+ ----------------------------------------------------------------------*/
+
+
+/* -------------------------------------------------------------------- */
+
+double PairSMTBQ::Primfcoup2(double c, double x, double delta)
+{
+
+ return (c*(c*c*c - 4* c*c* x - 4* (x - 2 *delta) * (x+delta)*(x+delta) +
+ 6* c *(x*x - delta*delta)))/(16* delta*delta*delta);
+
+}
+
+
+/* ----------------------------------------------------------------------
+ Integral of cutting function for derive (force) between x and c
+ ----------------------------------------------------------------------*/
+
+
+/* -------------------------------------------------------------------- */
+
+double PairSMTBQ::Intfcoup2(double c, double x, double delta)
+{
+
+ if (c x+delta)
+ {return Primfcoup2(x+delta,x,delta) - Primfcoup2(x,x,delta) ;}
+ else
+ {
+ return Primfcoup2(c,x,delta) - Primfcoup2(x,x,delta) ;}
+}
+
+
+/* ---------------------------------------------------------------------
+ Energy derivation respect charge Q
+ --------------------------------------------------------------------- */
+
+void PairSMTBQ::QForce_charge(int loop)
+{
+ int i,j,ii,jj,jnum;
+ int itype,jtype,m,gp;
+ double xtmp,ytmp,ztmp,evdwlCoul,fpairCoul;
+ double rsq,delr[3],chi_ij;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ double iq,jq,fqi,fqj,fqij,fqij2,fqjj;
+ int eflag;
+
+
+ double **x = atom->x;
+ double *q = atom->q;
+ int *type = atom->type;
+ int step = update->ntimestep;
+
+ int inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ int nlocal = atom->nlocal;
+
+
+ // loop over full neighbor list of my atoms
+
+ fqi = fqj = fqij = fqij2 = fqjj = 0.0;
+
+
+ // ==================
+ if (loop == 0) {
+ // ==================
+
+
+
+ for (ii = 0; ii < inum; ii ++) {
+ //--------------------------------
+ i = ilist[ii];
+ itype = map[type[i]];
+
+ gp = flag_QEq[i];
+
+ sbcov[i] =coord[i]= sbmet[i] = 0.0;
+
+ itype = map[type[i]];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ iq = q[i];
+
+
+
+
+ // two-body interactions
+
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+// -------------------------------
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ jtype = map[type[j]];
+ jq = q[j];
+
+ m = intype[itype][jtype];
+ if (intparams[m].intsm == 0) continue ;
+
+
+ delr[0] = x[j][0] - xtmp;
+ delr[1] = x[j][1] - ytmp;
+ delr[2] = x[j][2] - ztmp;
+ rsq = vec3_dot(delr,delr);
+
+
+// Covalente charge forces - sbcov initialization
+// ------------------------------------------------
+ if (sqrt(rsq) > intparams[m].dc2) continue;
+
+ attractive (&intparams[m],rsq,eflag,i,iq,j,jq);
+
+
+ } // ---------------------------------------------- for jj
+
+
+ } // -------------------------------------------- ii
+
+ // ===============================================
+ // Communicates the tables *sbcov and *sbmet
+ // to calculate the N-body forces
+ // ================================================
+
+ forward (sbcov) ; reverse (sbcov);
+ forward (coord) ; reverse (coord);
+ forward (sbmet) ; reverse (sbmet);
+
+
+ if (nteam == 0) return; //no oxide
+ if (Qstep == 0 || fmod(double(step), double(Qstep)) != 0.0) return;
+
+ // =======================
+ } // End of If(loop)
+ // =======================
+
+
+ // ===============================================
+
+ for (ii=0; ii cutmax) continue;
+
+ // 1/r charge forces
+ // --------------------
+ fqij = 0.0;
+// pot_ES2 (i,j,rsq,fqij2);
+ pot_ES (i,j,rsq,fqij);
+
+ potmad[i] += jq*fqij ;
+
+
+ } // ------ jj
+
+ fqi = 0.0;
+ pot_COV (¶ms[itype],i,fqi) ;
+ potcov[i] = fqi ;
+
+
+ } // ------- ii
+
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSMTBQ::Charge()
+{
+ int i,ii,iloop,itype,gp,m;
+ int *ilist;
+ double heatpq,qmass,dtq,dtq2,qtot,qtotll;
+ double t_init,t_end,dt;
+
+ double *Transf,*TransfAll;
+
+ double *q = atom->q;
+ double **x = atom->x;
+ int *type = atom->type;
+ int step = update->ntimestep;
+
+ int inum = list->inum;
+ ilist = list->ilist;
+
+
+ if (me == 0) t_init = MPI_Wtime();
+ if (step == 0) cluster = 0;
+
+ // ---------------------------
+ // Mise en place des groupes
+ // ---------------------------
+
+ if (strcmp(QEqMode,"BulkFromSlab") == 0)
+ groupBulkFromSlab_QEq();
+ else if (strcmp(QEqMode,"QEqAll") == 0)
+ groupQEqAll_QEq();
+ else if (strcmp(QEqMode,"QEqAllParallel") == 0)
+ groupQEqAllParallel_QEq();
+ else if (strcmp(QEqMode,"Surface") == 0)
+ groupSurface_QEq();
+
+
+
+ if (nteam+1 != cluster) {
+ memory->destroy(nQEqall);
+ memory->destroy(nQEqaall);
+ memory->destroy(nQEqcall);
+
+ cluster = nteam+1;
+ memory->create(nQEqall,nteam+1,"pair:nQEq");
+ memory->create(nQEqaall,nteam+1,"pair:nQEqa");
+ memory->create(nQEqcall,nteam+1,"pair:nQEqc");
+ }
+ // ---------------------------
+
+
+ double enegtotall[nteam+1],enegchkall[nteam+1],enegmaxall[nteam+1],qtota[nteam+1],qtotc[nteam+1];
+ double qtotcll[nteam+1],qtotall[nteam+1];
+ double sigmaa[nteam+1],sigmac[nteam+1],sigmaall[nteam+1],sigmacll[nteam+1];
+ int end[nteam+1], nQEq[nteam+1],nQEqc[nteam+1],nQEqa[nteam+1];
+
+
+ iloop = 0;
+
+
+ heatpq = 0.07;
+ qmass = 0.000548580;
+ dtq = 0.0006; // 0.0006
+ dtq2 = 0.5*dtq*dtq/qmass;
+
+ double enegchk[nteam+1];
+ double enegtot[nteam+1];
+ double enegmax[nteam+1];
+
+
+
+ for (i=0; intimestep == 0 && (strcmp(QInitMode,"true") == 0) ) {
+ //Carefull here it won't be fine if there are more than 2 species!!!
+ QOxInit=max(QOxInit, -0.98* params[1].qform *nQEqcall[gp]/nQEqaall[gp]) ;
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii]; itype = map[type[i]];
+ gp = flag_QEq[i];
+
+// if (gp == 0) continue;
+
+ if (itype == 0) q[i] = QOxInit ;
+ if (itype > 0) q[i] = -QOxInit * nQEqaall[gp]/nQEqcall[gp];
+ }
+ }
+
+ if (nteam == 0 || Qstep == 0) return;
+ if (fmod(double(step), double(Qstep)) != 0.0) return;
+ // --------------------------------------
+
+ // ----
+
+ // ----
+
+
+ if (me == 0 && strcmp(Bavard,"false") != 0) {
+ for (gp = 0; gp < nteam+1; gp++) {
+ printf (" ++++ Groupe %d - Nox %d Ncat %d\n",gp,nQEqaall[gp],nQEqcall[gp]);
+ if (nQEqcall[gp] !=0 && nQEqaall[gp] !=0 )
+ printf (" neutralite des charges %f\n qtotc %f qtota %f\n",
+ qtotll,qtotcll[gp]/nQEqcall[gp],qtotall[gp]/nQEqaall[gp]);
+ printf (" ---------------------------- \n");}
+ }
+
+ // ======================= Tab transfert ==================
+ // Transf[gp] = enegtot[gp]
+ // Transf[gp+cluster] = Qtotc[gp]
+ // Transf[gp+2cluster] = Qtota[gp]
+ // Transf[3cluster] = Qtot
+ // -------------------------------------------------------
+ memory->create(Transf,3*cluster+1,"pair:Tranf");
+ memory->create(TransfAll,3*cluster+1,"pair:TranfAll");
+ // ========================================================
+
+
+ // --------------------------------------------
+ for (iloop = 0; iloop < loopmax; iloop ++ ) {
+ // --------------------------------------------
+
+ qtot = qtotll = Transf[3*cluster] = 0.0 ;
+ for (gp=0; gpdestroy(Transf);
+ memory->destroy(TransfAll);
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSMTBQ::groupBulkFromSlab_QEq()
+{ int ii,i;
+ double **x=atom->x;
+ int *ilist;
+ double ztmp;
+ int inum = list->inum;
+ ilist = list->ilist;
+
+ for (ii = 0; ii < inum; ii++)
+ {
+ i = ilist[ii];
+ ztmp = x[i][2];
+ if (ztmp>zlim1QEq && ztmp< zlim2QEq)
+ flag_QEq[i]=1;
+ else
+ flag_QEq[i]=0;
+
+ nteam=1;
+
+ }
+
+}
+
+// ----------------------------------------------
+
+void PairSMTBQ::groupQEqAll_QEq()
+{ int ii,i;
+ int *ilist;
+ int inum = list->inum;
+ ilist = list->ilist;
+
+ nteam=1;
+
+ for (ii = 0; ii < inum; ii++)
+ {
+ i= ilist[ii];
+ flag_QEq[i]=1;
+ }
+
+}
+
+// ----------------------------------------------
+
+void PairSMTBQ::groupSurface_QEq()
+{ int ii,i;
+ double **x=atom->x;
+ int *ilist;
+ double ztmp;
+ int inum = list->inum;
+ ilist = list->ilist;
+
+ for (ii = 0; ii < inum; ii++)
+ {
+ i = ilist[ii];
+ ztmp = x[i][2];
+ if (ztmp>zlim1QEq)
+ flag_QEq[i]=1;
+ else
+ flag_QEq[i]=0;
+
+ nteam=1;
+
+ }
+
+}
+
+
+void PairSMTBQ::groupQEqAllParallel_QEq()
+{
+ int ii,i,jj,j,kk,k,itype,jtype,ktype,jnum,m,gp,zz,z,kgp;
+ int iproc,team_elt[10][nproc],team_QEq[10][nproc][5];
+ int *ilist,*jlist,*numneigh,**firstneigh,ngp,igp,nboite;
+ double delr[3],xtmp,ytmp,ztmp,rsq;
+ int **flag_gp, *nelt, **tab_gp;
+ int QEq,QEqall[nproc];
+
+ double **x = atom->x;
+ int *type = atom->type;
+ int inum = list->inum;
+ int nlocal = atom->nlocal;
+ int nghost = atom->nghost;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+
+ // +++++++++++++++++++++++++++++++++++++++++++++++++
+ // On declare et initialise nos p'tits tableaux
+ // +++++++++++++++++++++++++++++++++++++++++++++++++
+
+ nboite = nlocal + nghost;
+ int **tabtemp,**Alltabtemp, *gptmp, *Allgptmp;
+
+ memory->create(tabtemp,10*nproc+10,nproc,"pair:tabtemp");
+ memory->create(Alltabtemp,10*nproc+10,nproc,"pair:Alltabtemp");
+ memory->create(gptmp,10*nproc+10,"pair:gptmp");
+ memory->create(Allgptmp,10*nproc+10,"pair:Allgptmp");
+
+ memory->create(flag_gp,nproc,nboite,"pair:flag_gp");
+ memory->create(nelt,nboite,"pair:nelt");
+ memory->create(tab_gp,10,nboite,"pair:flag_gp");
+
+
+ for (i = 0; i < nlocal+nghost ; i++) { flag_QEq[i] = 0; }
+ for (i = 0; i < 10*nproc; i++) {
+ gptmp[i] = 0; Allgptmp[i] = 0;
+ for (j=0;jdestroy(tabtemp);
+ memory->destroy(Alltabtemp);
+ memory->destroy(gptmp);
+ memory->destroy(Allgptmp);
+ memory->destroy(flag_gp);
+ memory->destroy(tab_gp);
+ memory->destroy(nelt);
+
+ return;
+ }
+ // :::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+
+ for (m = 0; m < nproc; m++) {
+ for (i = 0; i < nboite; i++) { flag_gp[m][i] = 0; }
+ }
+
+ // OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
+ // It includes oxygens entering the QEq scheme O
+ // OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
+
+
+ ngp = igp = 0; nelt[ngp] = 0;
+
+ // On prend un oxygène
+ // printf ("[me %d] On prend un oxygene\n",me);
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii] ; itype = map[type[i]];
+ if (itype != 0 || flag_QEq[i] == 0) continue;
+
+ m = 0;
+
+ if (ngp != 0 && flag_gp[me][i] == ngp) continue;
+
+
+ // Grouping Initialisation
+ // -----------------------------
+ if (flag_gp[me][i] == 0) {
+ ngp += 1; nelt[ngp] = 0;
+ tab_gp[ngp][nelt[ngp]] = i;
+ flag_gp[me][i] = ngp;
+ nelt[ngp] += 1;
+ }
+ // -------------------------------
+
+
+ // Loop on the groups
+ // ----------------------
+ for (kk = 0; kk < nelt[ngp]; kk++)
+ {
+ k = tab_gp[ngp][kk];
+ ktype = map[type[k]];
+ // printf ("[me %d] kk - gp %d elemt %d : atom %d(%d)\n",me,ngp,kk,k,ktype);
+ if (k >= nlocal) continue;
+
+ xtmp = x[k][0];
+ ytmp = x[k][1];
+ ztmp = x[k][2];
+
+ // Loop on the oxygen's neighbor of the group
+ // ---------------------------------------------
+ jlist = firstneigh[k];
+ jnum = numneigh[k];
+ for (j = 0; j < nboite; j++ )
+ {
+ jtype = map[type[j]];
+ if (jtype == ktype) continue;
+ m = intype[itype][jtype];
+
+ if (jtype == 0 && flag_QEq[j] == 0) continue;
+
+
+ if (flag_gp[me][j] == ngp) continue;
+
+ delr[0] = x[j][0] - xtmp;
+ delr[1] = x[j][1] - ytmp;
+ delr[2] = x[j][2] - ztmp;
+ rsq = vec3_dot(delr,delr);
+
+ // -------------------------------------
+ if (sqrt(rsq) <= cutmax) {
+
+ flag_QEq[j] = 1; //Entre dans le schema QEq
+
+ // :::::::::::::::::::: Meeting of two group in the same proc :::::::::::::::::::::
+
+ if (flag_gp[me][j] != 0 && flag_gp[me][j] != ngp && nelt[flag_gp[me][j]] != 0) {
+ printf("[me %d] (atom %d) %d [elt %d] rencontre un nouveau groupe %d [elt %d] (atom %d)\n",
+ me,k,ngp,nelt[ngp],flag_gp[me][j],nelt[flag_gp[me][j]],j);
+
+ // On met a jours les tableaux
+ // -----------------------------
+ igp = flag_gp[me][j];
+ z = min(igp,ngp);
+
+ if (z == igp) { igp = z; }
+ else if (z == ngp) {
+ ngp = igp ; igp = z;
+ flag_gp[me][j] = ngp;
+ }
+
+ for (zz = 0; zz < nelt[ngp]; zz++) {
+ z = tab_gp[ngp][zz];
+ tab_gp[igp][nelt[igp]] = z;
+ nelt[igp] += 1;
+ flag_gp[me][z] = igp;
+ tab_gp[ngp][zz] = 0;
+ }
+
+ nelt[ngp] = 0;
+ for (z = nlocal; z < nboite; z++) {
+ if (flag_gp[me][z] == ngp) flag_gp[me][z] = igp;
+ }
+
+ m = 1; kk = 0;
+ ngp = igp;
+ break;
+ }
+ // ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+ flag_gp[me][j] = ngp;
+ if (j < nlocal)
+ {
+ tab_gp[ngp][nelt[ngp]] = j;
+ nelt[ngp] += 1;
+ }
+ }
+ } // for j
+ } // for k
+ } // for ii
+
+ // OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
+ // Groups communication
+ // OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
+ for (i = 0; i < nproc; i++) {
+ forward_int(flag_gp[i]); reverse_int(flag_gp[i]);
+ }
+ // ---
+
+
+ // =======================================================
+ // Loop on the cation to make them joined in the oxygen's
+ // group which it interacts
+ // =======================================================
+ igp = 0;
+ for (ii = 0; ii < inum; ii++)
+ {
+ i = ilist[ii] ; itype = map[type[i]];
+ if (itype == 0) continue;
+
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+ for (jj = 0; jj < jnum; jj++ )
+ {
+ j = jlist[jj] ; jtype = map[type[j]];
+ if (jtype != 0) continue;
+
+ m = 0;
+ for (iproc = 0; iproc < nproc; iproc++) {
+ if (flag_gp[iproc][j] != 0) m = flag_gp[iproc][j];
+ }
+ if (m == 0) continue;
+
+ delr[0] = x[j][0] - xtmp;
+ delr[1] = x[j][1] - ytmp;
+ delr[2] = x[j][2] - ztmp;
+ rsq = vec3_dot(delr,delr);
+
+ // ----------------------------------------
+ if (sqrt(rsq) <= cutmax) {
+ // if (sqrt(rsq) <= intparams[m].dc2) {
+ // ----------------------------------------
+
+ flag_QEq[i] = 1; igp = flag_gp[me][j];
+
+ if (flag_gp[me][i] == 0) flag_gp[me][i] = igp;
+
+ if (flag_gp[me][i] != igp && igp != 0) {
+ printf ("[me %d] Cation i %d gp %d [nelt %d] rencontre j %d(%d)du groupe %d [nelt %d]\n",
+ me,i,flag_gp[me][i],nelt[flag_gp[me][i]],j,jtype,igp,nelt[igp]);
+
+ igp = min(flag_gp[me][i],flag_gp[me][j]);
+ if (igp == flag_gp[me][i]) { kgp = flag_gp[me][j]; }
+ else { kgp = flag_gp[me][i]; }
+
+ for (k = 0; k < nelt[kgp]; k++) {
+ z = tab_gp[kgp][k];
+ tab_gp[igp][nelt[igp]] = z;
+ nelt[igp] += 1;
+ flag_gp[me][z] = igp;
+ tab_gp[kgp][k] = 0;
+ }
+ nelt[kgp] = 0;
+
+ for (k = 0; k < nboite; k++) {
+ if (flag_gp[me][k] == kgp) flag_gp[me][k] = igp;
+ }
+
+ }
+ m = 0;
+ for (k = 0; k < nelt[igp]; k++) {
+ if (tab_gp[igp][k] == i) m = 1;
+ }
+
+ if (i >= nlocal || m == 1 ) continue;
+ // printf ("[me %d] igp %d - nelt %d atom %d\n",me,igp,nelt[igp],i);
+ tab_gp[igp][nelt[igp]] = i;
+ nelt[igp] += 1;
+ break;
+ }
+
+ } // voisin j
+
+ } // atom i
+
+ /* ==================================================
+ Group Communication between proc for unification
+ ================================================== */
+ for (i = 0; i < nproc; i++) {
+ forward_int(flag_gp[i]); reverse_int(flag_gp[i]);
+ }
+
+ // =============== End of COMM =================
+
+
+ for (i = 0; i < nboite; i++) {
+
+ m = 10*me + flag_gp[me][i];
+ if (m == 10*me) continue; // Pas de groupe zero
+ gptmp[m] = 1;
+ for (k = 0; k < nproc; k++) {
+
+ if (k == me) continue;
+ if (tabtemp[m][k] != 0) continue;
+
+ if (flag_gp[k][i] != 0) {
+ tabtemp[m][k] = 10*k + flag_gp[k][i];
+ }
+ }
+
+ }
+
+ for (k = 0; k < 10*nproc; k++) {
+ MPI_Allreduce(tabtemp[k],Alltabtemp[k],nproc,MPI_INT,MPI_SUM,world); }
+ MPI_Allreduce(gptmp,Allgptmp,10*nproc,MPI_INT,MPI_SUM,world);
+
+ nteam = 0; iproc = 0;
+ for (igp = 0; igp < 10*nproc; igp++) {
+ if (Allgptmp[igp] == 0) continue;
+ iproc = int(double(igp)/10.0);
+ ngp = igp - 10*iproc;
+ if (nteam == 0) {
+
+ nteam += 1;
+ team_elt[nteam][iproc] = 0;
+ team_QEq[nteam][iproc][team_elt[nteam][iproc]] = ngp;
+ team_elt[nteam][iproc] += 1;
+ } else {
+ m = 0;
+ for (i = 1; i < nteam+1; i++) {
+ for (k = 0; k < team_elt[i][iproc]; k++) {
+ if (ngp == team_QEq[i][iproc][k]) m = 1;
+ } }
+ if (m == 1) continue;
+ // create a new team!!
+ // ---------------------------
+ if (m == 0) {
+ nteam += 1;
+ team_elt[nteam][iproc] = 0;
+ team_QEq[nteam][iproc][team_elt[nteam][iproc]] = ngp;
+ team_elt[nteam][iproc] += 1;
+ }
+ }
+ // -------
+ // On a mis une graine dans le groupe et nous allons
+ // remplir se groupe en questionnant le "tabtemp[igp][iproc]=jgp"
+ //
+ for (kk = 0; kk < nproc; kk++) {
+ for (k = 0; k < team_elt[nteam][kk]; k++) {
+
+ // On prend le gp le k ieme element de la team nteam sur le proc iproc
+ // ngp = 0;
+ ngp = team_QEq[nteam][kk][k];
+ kgp = 10*kk + ngp;
+
+ // On regarde sur les autre proc si ce gp ne pointe pas vers un autre.
+ for (i = 0; i < nproc; i++) {
+ if (i == kk) continue;
+ if (Alltabtemp[kgp][i] == 0) continue;
+
+ if (Alltabtemp[kgp][i] != 0) ngp = Alltabtemp[kgp][i];
+ ngp = ngp - 10*i;
+
+ // Est ce que ce groupe est nouveau?
+ m = 0;
+ for (j = 0; j < team_elt[nteam][i]; j++) {
+ if (team_QEq[nteam][i][j] == ngp) m = 1;
+ }
+
+ if (m == 0) {
+ iproc = i; k = 0;
+ team_QEq[nteam][i][team_elt[nteam][i]] = ngp ;
+ team_elt[nteam][i] += 1;
+ }
+ } // regard sur les autre proc
+
+ } // On rempli de proche en proche
+ } // boucle kk sur les proc
+ }
+
+ // Finalement on met le numero de la team en indice du flag_QEq, c mieu!
+ // ---------------------------------------------------------------------
+
+ for (ii = 0; ii < inum; ii++)
+ {
+ i = ilist[ii]; m = 0; itype = map[type[i]];
+ if (flag_QEq[i] == 0) continue;
+
+ gp = flag_gp[me][i];
+ for (j = 1; j <= nteam; j++) {
+ for (k = 0; k < team_elt[j][me]; k++) {
+ if (gp == team_QEq[j][me][k]) {
+ flag_QEq[i] = j; m = 1;
+ break;
+ }
+ }
+ if (m == 1) break;
+ }
+ }
+
+ memory->destroy(tabtemp);
+ memory->destroy(Alltabtemp);
+ memory->destroy(gptmp);
+ memory->destroy(Allgptmp);
+ memory->destroy(flag_gp);
+ memory->destroy(tab_gp);
+ memory->destroy(nelt);
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSMTBQ::Init_charge(int *nQEq, int *nQEqa, int *nQEqc)
+{
+ int ii,i,gp,itype;
+ int *ilist,test[nteam],init[nteam];
+ double bound,tot,totll;
+
+ int inum = list->inum;
+ int *type = atom->type;
+ double *q = atom->q;
+ ilist = list->ilist;
+
+ if (nteam == 0) return;
+
+ if (me == 0) printf (" ======== Init_charge ======== \n");
+ for (gp = 0; gp < cluster; gp++) {
+ test[gp] = 0; init[gp] = 0;
+ }
+
+
+ // On fait un test sur les charges pour voir sont
+ // elles sont dans le domaine delimiter par DeltaQ
+ // -------------------------------------------------
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii]; itype = map[type[i]];
+ gp = flag_QEq[i];
+
+ if (gp != 0 && itype == 0) {
+ bound = fabs(2.0*ncov/params[itype].sto - fabs(params[itype].qform)) ;
+ if (bound == fabs(params[itype].qform)) continue;
+ if (fabs(q[i]) < bound) test[gp] = 1;
+ }
+ }
+
+ MPI_Allreduce(test,init,nteam+1,MPI_INT,MPI_SUM,world);
+
+ // On fait que sur les atomes hybrides!!!
+ // ----------------------------------------
+
+ tot = totll = 0.0;
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii]; itype = map[type[i]];
+ gp = flag_QEq[i];
+
+ if (gp != 0 && init[gp] != 0) {
+ if (itype == 0) q[i] = -1.96;
+ if (itype != 0) q[i] = 1.96 * double(nQEqaall[gp]) / double(nQEqcall[gp]);
+ }
+ tot += q[i];
+ }
+ MPI_Allreduce(&tot,&totll,1,MPI_INT,MPI_SUM,world);
+ if (me == 0) printf (" === Fin de init_charge qtot %20.15f ====\n",totll);
+
+}
+/* ----------------------------------------------------------------------
+ * COMMUNICATION
+ * ---------------------------------------------------------------------- */
+
+int PairSMTBQ::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i ++) {
+ j = list[i];
+ buf[m++] = tab_comm[j];
+// if (j < 3) printf ("[%d] %d pfc %d %d buf_send = %f \n",me,n,i,m-1,buf[m-1]);
+ }
+ return 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSMTBQ::unpack_forward_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n ;
+ for (i = first; i < last; i++) {
+ tab_comm[i] = buf[m++];
+// if (inlocal;
+ int nghost = atom->nghost;
+
+ for (i=0; iforward_comm_pair(this);
+
+ for (i=0; inlocal;
+ int nghost = atom->nghost;
+
+ for (i=0; ireverse_comm_pair(this);
+
+ for (i=0; inlocal;
+ int nghost = atom->nghost;
+
+ for (i=0; iforward_comm_pair(this);
+
+ for (i=0; i 0.1) tab[i] = int(tab_comm[i]) ; }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSMTBQ::reverse_int(int *tab)
+{
+ int i;
+ int nlocal = atom->nlocal;
+ int nghost = atom->nghost;
+
+ for (i=0; ireverse_comm_pair(this);
+
+ for (i=0; i 0.1) tab[i] = int(tab_comm[i]); }
+}
+
+/* ---------------------------------------------------------------------- */
+/* ---------------------------------------------------------------------- */
+
+double PairSMTBQ::memory_usage()
+{
+ double bytes = maxeatom * sizeof(double);
+ bytes += maxvatom*6 * sizeof(double);
+ bytes += nmax * sizeof(int);
+ bytes += MAXNEIGH * nmax * sizeof(int);
+ return bytes;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSMTBQ::add_pages(int howmany)
+{
+ int toppage = maxpage;
+ maxpage += howmany*PGDELTA;
+
+ pages = (int **)
+ memory->srealloc(pages,maxpage*sizeof(int *),"pair:pages");
+ for (int i = toppage; i < maxpage; i++)
+ memory->create(pages[i],pgsize,"pair:pages[i]");
+}
+
+/* ---------------------------------------------------------------------- */
+
+
+int PairSMTBQ::Tokenize( char* s, char*** tok )
+{
+ char test[MAXLINE];
+ const char *sep = "' ";
+ char *mot;
+ int count=0;
+ mot = NULL;
+
+
+ strncpy( test, s, MAXLINE );
+
+ for( mot = strtok(test, sep); mot; mot = strtok(NULL, sep) ) {
+ strncpy( (*tok)[count], mot, MAXLINE );
+ count++;
+ }
+
+ return count;
+}
+
+
+
+void PairSMTBQ::CheckEnergyVSForce()
+{
+ double drL,iq,jq,rsq,evdwlCoul,fpairCoul,eflag,ErepR,frepR,fpair,evdwl;
+ int i,j,iiiMax,iii,iCoord;
+ int itype,jtype,l,m;
+ double r,t1,t2,sds,xi,engSurf,fforceSurf;
+ double eng,fforce,engBB,fforceBB;
+
+ double za,zb,gam,dgam,dza,dzb,
+ d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2,na,nb;
+ int *type = atom->type;
+ char *NameFile;
+
+
+ i=0;
+ j=1;
+ map[type[i]]=0; //ox
+ itype=map[type[i]];
+ iq=-1;
+
+
+ map[type[j]]=0; //ox
+ jtype=map[type[j]];
+ jq=-1;
+ coord[i]=coordOxBulk;
+ coord[j]=coordOxBulk;
+ m = intype[itype][jtype];
+
+
+ na = params[itype].ne ;
+ nb = params[jtype].ne ;
+ za = params[itype].dzeta ;
+ zb = params[jtype].dzeta ;
+
+
+ // Ouverture du fichier
+
+ for (iCoord=1;iCoord<5; iCoord++)
+ {
+
+ if (iCoord==1)
+ {coord[i]=2.2;
+ coord[j]=2.1;
+ NameFile="energyandForceOxOxUnderCoord.txt";
+ }
+ if (iCoord==2)
+ {coord[i]=coordOxBulk;
+ coord[j]=coordOxBulk;
+ NameFile="energyandForceOxOxCoordBulk.txt";
+ }
+ if (iCoord==3)
+ {coord[i]=3.8;
+ coord[j]=4;
+ NameFile="energyandForceOxOxOverCoord.txt";
+ }
+ if (iCoord==4)
+ {coord[i]=2.2;
+ coord[j]=3.5;
+ NameFile="energyandForceOxOxUnderOverCoord.txt";
+ }
+
+
+ ofstream fichierOxOx(NameFile, ios::out | ios::trunc) ;
+
+ drL=0.0001;
+ iiiMax=int((cutmax-1.2)/drL);
+ for (iii=1; iii< iiiMax ; iii++){
+ r=1.2+drL*iii;
+ rsq=r*r;
+ evdwlCoul = 0.0 ; fpairCoul = 0.0;
+ potqeq(i,j,iq,jq,rsq,fpairCoul,eflag,evdwlCoul);
+ fpairCoul=fpairCoul*r;
+
+ rep_OO (&intparams[m],rsq,fpair,eflag,evdwl);
+ ErepR = evdwl;
+ frepR= fpair*r;
+
+ gam = dgam = dza = dzb = d2zaa = d2zab =
+ d2zbb = d2zra = d2zrb = d2gamr2 = 0.0 ;
+
+
+// gammas_(na,nb,za,zb,r,gam,dgam,dza,dzb,
+// d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
+ gammas(na,nb,za,zb,r,gam,dgam,dza,dzb,
+ d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
+
+
+ sds = rsq/ds ; l = int(sds) ;
+ xi = sds - double(l) ;
+
+
+ t1 = fafb[l][m] + (fafb[l+1][m] - fafb[l][m])*xi;
+ t2 = fafb[l+1][m] + (fafb[l+2][m] - fafb[l+1][m])*(xi-1.0);
+ eng = iq*jq*(t1 + (t2 - t1)*xi/2.0);
+
+ t1 = dfafb[l][m] + (dfafb[l+1][m] - dfafb[l][m])*xi;
+ t2 = dfafb[l+1][m] + (dfafb[l+2][m] - dfafb[l+1][m])*(xi-1);
+ fforce = - iq*jq*(t1 + (t2 - t1)*xi/2.0)*r ;
+
+
+ t1 = fafbOxOxSurf[l] + (fafbOxOxSurf[l+1] - fafbOxOxSurf[l])*xi;
+ t2 = fafbOxOxSurf[l+1] + (fafbOxOxSurf[l+2] - fafbOxOxSurf[l+1])*(xi-1.0);
+ engSurf = iq*jq*(t1 + (t2 - t1)*xi/2.0);
+
+ t1 = fafbOxOxBB[l] + (fafbOxOxBB[l+1] - fafbOxOxBB[l])*xi;
+ t2 = fafbOxOxBB[l+1] + (fafbOxOxBB[l+2] - fafbOxOxBB[l+1])*(xi-1.0);
+ engBB = iq*jq*(t1 + (t2 - t1)*xi/2.0);
+
+ t1 = dfafbOxOxSurf[l] + (dfafbOxOxSurf[l+1] - dfafbOxOxSurf[l])*xi;
+ t2 = dfafbOxOxSurf[l+1] + (dfafbOxOxSurf[l+2] - dfafbOxOxSurf[l+1])*(xi-1);
+ fforceSurf = - iq*jq*(t1 + (t2 - t1)*xi/2.0)*r ;
+
+ t1 = dfafbOxOxBB[l] + (dfafbOxOxBB[l+1] - dfafbOxOxBB[l])*xi;
+ t2 = dfafbOxOxBB[l+1] + (dfafbOxOxBB[l+2] - dfafbOxOxBB[l+1])*(xi-1);
+ fforceBB = - iq*jq*(t1 + (t2 - t1)*xi/2.0)*r ;
+
+ if (fichierOxOx) { fichierOxOx<< setprecision (9) <= 1; i--) {
+ rtrm=rtrm*halfr;
+ drtrm=drtrm*halfr;
+ ztrm=ztrm*2.0*zb;
+ ctrm=ztrm/factorial(nb2-i);
+
+ css(ss,na2-1,nb2-i,z2ra,z2rb,r,deriv,dzeta1,dzeta2,
+ d2zeta11,d2zeta12,d2zeta22,d2zeta1r,d2zeta2r,deriv2);
+
+ trm1=float(i)*ctrm;
+ trm2=trm1*rtrm;
+ gam=gam-trm2*ss;
+ trm3=trm1*float(na2+nb2-i)*drtrm;
+ dgam=dgam-trm2*deriv-0.5*trm3*ss;
+ d2gamr2=d2gamr2-trm2*deriv2-trm3*deriv-0.5*trm3*float(na2+nb2-i-1)*ss/r;
+ dza=dza-trm2*dzeta1;
+ dzb=dzb-(trm2/zb)*((float(nb2-i))*ss+zb*dzeta2);
+ d2zaa=d2zaa-trm2*d2zeta11;
+ d2zab=d2zab-(trm2/zb)*((float(nb2-i))*dzeta1+zb*d2zeta12);
+ d2zbb=d2zbb-(trm2/zb)*(2.0*(float(nb2-i))*dzeta2+zb*d2zeta22 +
+ (float((nb2-i-1)*(nb2-i))*ss/zb));
+ d2zra=d2zra-trm2*d2zeta1r-0.5*trm3*dzeta1;
+ d2zrb=d2zrb-(trm2/zb)*((float(nb2-i))*deriv+zb*d2zeta2r) -
+ 0.5*(trm3/zb)*((float(nb2-i))*ss+zb*dzeta2);
+ }
+
+// Multiply by coefficients
+
+ trm3=pow(2.0*za,na2+1)/factorial(na2);
+ gam=gam*trm3;
+ dgam=dgam*trm3;
+ rfct1=((float(na2+1))/za);
+ rgam1=rfct1*gam;
+ dza=dza*trm3;
+ rdza1=2.0*rfct1*dza;
+ dza=dza+rgam1;
+ dzb=dzb*trm3;
+ rgam2=rgam1*float(na2)/za;
+ d2zaa=d2zaa*trm3+rgam2+rdza1;
+ d2zab=d2zab*trm3+rfct1*dzb;
+ d2zbb=d2zbb*trm3;
+ d2zra=d2zra*trm3+rfct1*dgam;
+ d2zrb=d2zrb*trm3;
+ d2gamr2=d2gamr2*trm3;
+ return;
+
+ }
+/* --------------------------------------------------------------------------------
+ Css
+ -------------------------------------------------------------------------------- */
+void PairSMTBQ::css(double &s, double nn1, double nn2, double alpha, double beta, double r,
+ double &deriv, double &dzeta1, double &dzeta2, double &d2zeta11, double &d2zeta12,
+ double &d2zeta22, double &d2zeta1r, double &d2zeta2r, double &deriv2)
+{
+// implicit real (a-h,o-z)
+// common /fctrl/ fct(30) // A RAJOUTER DANS Pair_SMTBQ.h
+/* ------------------------------------------------------------------
+ Modified integral calculation routine for Slater orbitals
+ including derivatives. This version is for S orbitals only.
+
+ dzeta1 and dzeta2 are the first derivatives with respect to zetas
+ and d2zeta11/d2zeta12/d2zeta22 are the second.
+ d2zeta1r and d2zeta2r are the mixed zeta/r second derivatives
+ deriv2 is the second derivative with respect to r
+
+ Julian Gale, Imperial College, December 1997
+ ------------------------------------------------------------------- */
+ int i,i1,nni1; //ulim
+ double ulim,n1,n2,p,pt,x,k,dpdz1,dpdz2,dptdz1,dptdz2,dpdr,dptdr,d2pdz1r,
+ d2pdz2r,d2ptdz1r,d2ptdz2r,zeta1,zeta2,sumzeta,difzeta,coff;
+
+ double da1[30],da2[30],db1[30],db2[30];
+ double d2a11[30],d2a12[30],d2a22[30],dar[30];
+ double d2b11[30],d2b12[30],d2b22[30],dbr[30];
+ double d2a1r[30],d2a2r[30],d2b1r[30],d2b2r[30];
+ double d2ar2[30],d2br2[30];
+
+ double *a,*b;
+ memory->create(a,31,"pair:a");
+ memory->create(b,31,"pair:a");
+
+
+// Set up factorials - stored as factorial(n) in location(n+1)
+
+ for (i = 1; i <= 30; i++) {
+ fct[i]=factorial(i-1);
+ }
+ dzeta1=0.0;
+ dzeta2=0.0;
+ d2zeta11=0.0;
+ d2zeta12=0.0;
+ d2zeta22=0.0;
+ d2zeta1r=0.0;
+ d2zeta2r=0.0;
+ deriv=0.0;
+ deriv2=0.0;
+ n1=nn1;
+ n2=nn2;
+ p =(alpha + beta)*0.5;
+ pt=(alpha - beta)*0.5;
+ x = 0.0;
+ zeta1=alpha/r;
+ zeta2=beta/r;
+ sumzeta=zeta1+zeta2;
+ difzeta=zeta1-zeta2;
+
+// Partial derivative terms for zeta derivatives
+
+ dpdz1=r;
+ dpdz2=r;
+ dptdz1=r;
+ dptdz2=-r;
+ dpdr=0.5*sumzeta;
+ dptdr=0.5*difzeta;
+ d2pdz1r=1.0;
+ d2pdz2r=1.0;
+ d2ptdz1r=1.0;
+ d2ptdz2r=-1.0;
+
+// Reverse quantum numbers if necessary -
+// also change the sign of difzeta to match
+// change in sign of pt
+
+ if (n2 < n1) {
+ k = n1;
+ n1= n2;
+ n2= k;
+ pt=-pt;
+ difzeta=-difzeta;
+ dptdr=-dptdr;
+ dptdz1=-dptdz1;
+ dptdz2=-dptdz2;
+ d2ptdz1r=-d2ptdz1r;
+ d2ptdz2r=-d2ptdz2r;
+ }
+
+// Trap for enormously long distances which would cause
+// caintgs or cbintgs to crash with an overflow
+
+ if (p > 86.0 || pt > 86.0) {
+ s=0.0;
+ return;
+ }
+//***************************
+// Find a and b integrals *
+//***************************
+ caintgs(p,n1+n2+3,a);
+ cbintgs(pt,n1+n2+3,b);
+
+// Convert derivatives with respect to p and pt
+// into derivatives with respect to zeta1 and
+// zeta2
+
+ ulim=n1+n2+1;
+ for (i = 1; i <= int(ulim); i++) {
+ da1[i]=-a[i+1]*dpdz1;
+ da2[i]=-a[i+1]*dpdz2;
+ db1[i]=-b[i+1]*dptdz1;
+ db2[i]=-b[i+1]*dptdz2;
+ d2a11[i]=a[i+2]*dpdz1*dpdz1;
+ d2a12[i]=a[i+2]*dpdz1*dpdz2;
+ d2a22[i]=a[i+2]*dpdz2*dpdz2;
+ d2b11[i]=b[i+2]*dptdz1*dptdz1;
+ d2b12[i]=b[i+2]*dptdz1*dptdz2;
+ d2b22[i]=b[i+2]*dptdz2*dptdz2;
+ dar[i]=-a[i+1]*dpdr;
+ dbr[i]=-b[i+1]*dptdr;
+ d2a1r[i]=a[i+2]*dpdz1*dpdr-a[i+1]*d2pdz1r;
+ d2a2r[i]=a[i+2]*dpdz2*dpdr-a[i+1]*d2pdz2r;
+ d2b1r[i]=b[i+2]*dptdz1*dptdr-b[i+1]*d2ptdz1r;
+ d2b2r[i]=b[i+2]*dptdz2*dptdr-b[i+1]*d2ptdz2r;
+ d2ar2[i]=a[i+2]*dpdr*dpdr;
+ d2br2[i]=b[i+2]*dptdr*dptdr;
+ }
+
+// Begin section used for overlap integrals involving s functions
+
+ for (i1 = 1; i1 <= int(ulim); i1++) {
+ nni1=n1+n2-i1+2;
+ coff=coeffs(n1,n2,i1-1);
+ deriv=deriv+coff*(dar[i1]*b[nni1]+a[i1]*dbr[nni1]);
+ x=x+coff*a[i1]*b[nni1];
+ dzeta1=dzeta1+coff*(da1[i1]*b[nni1]+a[i1]*db1[nni1]);
+ dzeta2=dzeta2+coff*(da2[i1]*b[nni1]+a[i1]*db2[nni1]);
+ d2zeta11=d2zeta11+coff*(d2a11[i1]*b[nni1]+a[i1]*d2b11[nni1]+
+ 2.0*da1[i1]*db1[nni1]);
+ d2zeta12=d2zeta12+coff*(d2a12[i1]*b[nni1]+a[i1]*d2b12[nni1]+
+ da1[i1]*db2[nni1]+da2[i1]*db1[nni1]);
+ d2zeta22=d2zeta22+coff*(d2a22[i1]*b[nni1]+a[i1]*d2b22[nni1]+
+ 2.0*da2[i1]*db2[nni1]);
+ d2zeta1r=d2zeta1r+coff*(d2a1r[i1]*b[nni1]+dar[i1]*db1[nni1]+
+ da1[i1]*dbr[nni1]+a[i1]*d2b1r[nni1]);
+ d2zeta2r=d2zeta2r+coff*(d2a2r[i1]*b[nni1]+dar[i1]*db2[nni1]+
+ da2[i1]*dbr[nni1]+a[i1]*d2b2r[nni1]);
+ deriv2=deriv2+coff*(d2ar2[i1]*b[nni1]+a[i1]*d2br2[nni1]+
+ 2.0*dar[i1]*dbr[nni1]);
+ }
+ s=x*0.5;
+ deriv=0.5*deriv;
+ deriv2=0.5*deriv2;
+ dzeta1=0.5*dzeta1;
+ dzeta2=0.5*dzeta2;
+ d2zeta11=0.5*d2zeta11;
+ d2zeta12=0.5*d2zeta12;
+ d2zeta22=0.5*d2zeta22;
+ d2zeta1r=0.5*d2zeta1r;
+ d2zeta2r=0.5*d2zeta2r;
+
+ memory->destroy(a);
+ memory->destroy(b);
+
+ return;
+}
+/* -------------------------------------------------------------------------------
+ coeffs
+ ------------------------------------------------------------------------------- */
+ double PairSMTBQ::coeffs(int na, int nb, int k)
+ {
+ // implicit real (a-h,o-z)
+ // common /fctrl/ fct(30)
+
+ int il,lie,je,ia,i,j,ie,l;
+ double coeffs;
+
+// Statement function
+// binm(n,i)=fct(n+1)/(fct(n-i+1)*fct(i+1));
+
+ coeffs=0.0;
+ l=na+nb-k;
+ ie=min(l,na)+1;
+ je=min(l,nb);
+ ia=l-je+1;
+ for (il = ia; il <= ie; il++) {
+ i=il-1;
+ j=l-i; // D'ou vient le i
+ coeffs=coeffs + binm(na,i)*binm(nb,j)*pow(-1.,j);
+ }
+ return coeffs;
+ }
+
+// ============================================
+
+double PairSMTBQ::binm(int n, int i)
+{
+ return fct[n+1]/(fct[n-i+1]*fct[i+1]);
+}
+
+/* ---------------------------------------------------------------------------------
+ Caintgs
+ -------------------------------------------------------------------------------- */
+ void PairSMTBQ::caintgs (double x, int k, double *a)
+ {
+// implicit real (a-h,o-z)
+// dimension a(30)
+ int i;
+ double cste,rx;
+
+ cste=exp(-x);
+ rx=1.0/x;
+ a[1]=cste*rx;
+ for (i = 1; i <= k; i++) {
+ a[i+1]=(a[i]*float(i)+cste)*rx;
+ }
+ return;
+ }
+/* -----------------------------------------------------------------------------------
+ Cbintgs
+ ----------------------------------------------------------------------------------- */
+void PairSMTBQ::cbintgs( double x, int k, double *b)
+{
+// implicit real (a-h,o-z)
+/* *******************************************************************
+! Fills array of b-integrals. note that b(i) is b(i-1) in the
+! usual notation
+! for x.gt.3 exponential formula is used
+! for 2.lt.x.le.3 and k.le.10 exponential formula is used
+! for 2.lt.x.le.3 and k.gt.10 15 term series is used
+! for 1.lt.x .e.2 and k.le.7 exponential formula is used
+! for 1.lt.x.le.2 and k.gt.7 12 term series is used
+! for .5.lt.x.le.1 and k.le.5 exponential formula is used
+! for .5.lt.x.le.1 and k.gt.5 7 term series is used
+! for x.le..5 6 term series is used
+!******************************************************************* */
+// dimension b(30)
+// common /fctrl/ fct(30)
+
+ int i0,m,last,i;
+ double absx,expx,expmx,ytrm,y,rx;
+
+ i0=0;
+ absx=fabs(x);
+
+ if (absx > 3.0) goto g120;
+ if (absx > 2.0) goto g20;
+ if (absx > 1.0) goto g50;
+ if (absx > 0.5) goto g80;
+ if (absx > 1.0e-8) goto g110;
+ goto g170;
+ g110: last=6;
+ goto g140;
+ g80: if (k <= 5) goto g120;
+ last=7;
+ goto g140;
+ g50: if (k <= 7) goto g120;
+ last=12;
+ goto g140;
+ g20: if (k <= 10) goto g120;
+ last=15;
+ goto g140;
+
+ g120: expx=exp(x);
+ expmx=1./expx;
+ rx=1.0/x;
+ b[1]=(expx-expmx)*rx;
+ for (i = 1; i <= k ; i++) {
+ b[i+1]=(float(i)*b[i]+ pow(-1.0,i)*expx-expmx)*rx;
+ }
+ goto g190;
+//
+// Series to calculate b(i)
+//
+ g140: for (i = i0; i <= k ; i++) {
+ y=0.;
+ for (m = i0; m <= last; m++) {
+ ytrm = pow(-x,m-1)*(1. - pow(-1.,m+i+1))/(fct[m+1]*float(m+i+1));
+ y = y + ytrm*(-x);
+ }
+ b[i+1] = y;
+ }
+ goto g190;
+//
+// x extremely small
+//
+ g170: for (i = i0; i <= k; i++) {
+ b[i+1] = (1.-pow(-1.,i+1))/float(i+1);
+ }
+ g190:
+ return;
+ }
+/* ------------------------------------------------------------------------------------
+ Factorial
+ ------------------------------------------------------------------------------------ */
+double PairSMTBQ::factorial(double n)
+{
+// implicit real(a-h,o-z)
+ double rn,factorial;
+ int i;
+
+//
+// Calculates the factorial of an integer n
+//
+ rn = n;
+ if (n <= 1.) {
+ factorial=1.0;
+ return factorial;
+ }
+ factorial=1.0;
+ for (i = 2; i <= int(n); i++) {
+ factorial = factorial*float(i);
+ }
+ return factorial;
+}
+
+
+/* ============================== This is the END... ================================== */
diff --git a/src/USER-SMTBQ/pair_smtbq.h b/src/USER-SMTBQ/pair_smtbq.h
new file mode 100644
index 0000000000..2a1dc3bb8a
--- /dev/null
+++ b/src/USER-SMTBQ/pair_smtbq.h
@@ -0,0 +1,197 @@
+/* ----------------------------------------------------------------------
+ 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(smtbq,PairSMTBQ)
+
+#else
+
+#ifndef LMP_PAIR_SMTBQ_H
+#define LMP_PAIR_SMTBQ_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairSMTBQ : public Pair {
+ public:
+ PairSMTBQ(class LAMMPS *);
+ virtual ~PairSMTBQ();
+ virtual void compute(int, int);
+
+ void settings(int, char **);
+ void coeff(int, char **);
+ void init_style();
+ double init_one(int, int);
+ double memory_usage();
+
+
+protected:
+ struct Param {
+ double sto,n0,ne,chi,dj;
+ double R,dzeta; //Rayon slater
+ double cutsq;
+ double qform, masse; // Charge formelle
+ char *nom;
+ };
+
+ struct Intparam {
+ double a,p,ksi,q,r0,dc1,dc2;
+ double abuck,rhobuck,neig_cut,aOO,bOO,r1OO,r2OO;
+ char *typepot,*mode;
+ int intsm;
+ };
+
+ double rmin,dr,ds; // table parameter
+ int kmax;
+ int Qstep; // Frenquence of charge resolution
+ double precision; // acurate of convergence
+ int loopmax; // max of interation
+
+ double cutmax; // max cutoff for all elements
+ int nelements; // # of unique elements
+ char **elements; // names of unique elements
+ char *QEqMode; // name of QEqMode
+ char *Bavard; // Verbose parameter
+ char *writepot; // write or not the electronegativity component
+ char *writeenerg; // write or not the energy component
+ char *QInitMode; // mode of initialization of charges
+ double zlim1QEq; // z limit for QEq equilibration
+ double zlim2QEq; // z limit for QEq equilibration
+ double QOxInit; // Initial charge for oxygen atoms (if the charge is not specified)
+ int *map; // mapping from atom types to elements
+ int nparams; // # of stored parameter sets
+ int maxparam; // max # of parameter sets
+ int maxintparam; // max # of interaction sets
+ int maxintsm; // max # of interaction SM
+ double r1Coord,r2Coord;
+ Param *params; // parameter set for an I atom
+ Intparam *intparams; // parameter set for an I interaction
+
+ int nmax,*nQEqall,*nQEqaall,*nQEqcall;
+ double *qf,*q1,*q2,Nevery,Neverypot;
+
+// Coulombian interaction
+ double *esm, **fafb, **dfafb, *fafbOxOxSurf, *dfafbOxOxSurf, *fafbTiOxSurf,*dfafbTiOxSurf;
+ double *potqn, *dpotqn, Vself, *Zsm,*coord, *fafbOxOxBB, *dfafbOxOxBB,*fafbTiOxBB, *dfafbTiOxBB ;
+ int **intype, **coultype;
+ int *NCo;
+ double coordOxBulk,coordOxSurf,ROxSurf,coordOxBB,ROxBB;
+
+// Covalent interaction
+ double *ecov, *potmad, *potself, *potcov, *chimet;
+ double **tabsmb,**dtabsmb, **tabsmr, **dtabsmr, *sbcov, *sbmet;
+ double ncov;
+
+// Neighbor Table
+ int nteam,cluster,*hybrid;
+ int *nvsm, **vsm, *flag_QEq;
+
+// Parallelisation
+ int me, nproc;
+ double *tab_comm;
+
+// GAMMAS function
+ double *fct;
+
+// HERE its routines
+// =====================================
+ void allocate();
+ virtual void read_file(char *);
+
+ void tabsm();
+ void tabqeq();
+
+ void potqeq(int, int, double, double, double,
+ double &, int, double &);
+ void pot_ES (int, int, double, double &);
+ void pot_ES2 (int, int, double, double &);
+
+ double self(Param *, double);
+ double qfo_self(Param *, double);
+
+ virtual void repulsive(Intparam *, double, int, int, double &, int, double &);
+ virtual void rep_OO (Intparam *, double, double &, int, double &);
+ virtual void Attr_OO (Intparam *, double, double &, int, double &);
+ virtual void attractive(Intparam *, double, int, int, double, int, double );
+ void f_att(Intparam *, int, int, double, double &) ;
+ void pot_COV(Param *, int, double &);
+ double potmet(Intparam *, double, int, double, int, double);
+
+ double fcoupure(double, double, double);
+ double fcoupured(double, double, double);
+
+ double fcoup2(double, double, double);
+ double Intfcoup2(double, double, double);
+ double Primfcoup2(double, double, double);
+
+ void groupBulkFromSlab_QEq();
+ void groupQEqAllParallel_QEq();
+ void groupQEqAll_QEq();
+ void groupSurface_QEq();
+ void QForce_charge(int);
+ void Charge();
+ void Init_charge(int*, int*, int*);
+ void CheckEnergyVSForce();
+
+// ===========================================
+// Communication pack
+ int pack_forward_comm (int, int*, double*, int, int*);
+ void unpack_forward_comm (int, int, double*);
+ int pack_reverse_comm (int, int, double*);
+ void unpack_reverse_comm (int, int*, double*);
+ void forward (double*); void reverse (double*);
+ void forward_int (int*); void reverse_int (int*);
+
+ int Tokenize( char* , char*** );
+
+ inline double vec3_dot(const double x[3], const double y[3]) const {
+ return x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
+ }
+
+ template const T& min ( const T& a, const T& b ) {
+ return !(b