From 0eb2e8f9cb98b70bd8aee0e4ef6d7ff16e2139f0 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 22 Oct 2015 22:56:21 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14172 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/Makefile | 2 +- src/USER-SMTBQ/README | 13 + src/USER-SMTBQ/pair_smtbq.cpp | 4241 +++++++++++++++++++++++++++++++++ src/USER-SMTBQ/pair_smtbq.h | 197 ++ 4 files changed, 4452 insertions(+), 1 deletion(-) create mode 100644 src/USER-SMTBQ/README create mode 100644 src/USER-SMTBQ/pair_smtbq.cpp create mode 100644 src/USER-SMTBQ/pair_smtbq.h 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