diff --git a/src/USER-SMTBQ/pair_smtbq.cpp b/src/USER-SMTBQ/pair_smtbq.cpp index 0955f7eee2..b2d17109be 100644 --- a/src/USER-SMTBQ/pair_smtbq.cpp +++ b/src/USER-SMTBQ/pair_smtbq.cpp @@ -9,7 +9,7 @@ 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 @@ -36,12 +36,12 @@ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details: . - ------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ -#include -#include -#include -#include +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" #include "pair_smtbq.h" #include "atom.h" #include "comm.h" @@ -52,32 +52,36 @@ #include "group.h" #include "update.h" #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "error.h" #include "domain.h" -#include -#include -#include #include +#include -//to DELETE -//#include -//#include -//#include -//#include -//to DELETE using namespace std; using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; #define MAXLINE 2048 #define MAXTOKENS 2048 #define DELTA 4 #define PGDELTA 1 #define MAXNEIGH 24 -#define Pi 4*atan(1) + +/* ------------------------------------------------------------------------------------ + + Calculates the factorial of an integer n via recursion. + + ------------------------------------------------------------------------------------ */ +static double factorial(int n) +{ + if (n <= 1) return 1.0; + else return static_cast(n)*factorial(n-1); +} /* ---------------------------------------------------------------------- */ @@ -282,13 +286,13 @@ void PairSMTBQ::coeff(int narg, char **arg) } 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++; - } + 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 @@ -374,7 +378,7 @@ double PairSMTBQ::init_one(int i, int j) void PairSMTBQ::read_file(char *file) { - int c, num_atom_types,i,k,m,test,j,verbose; + int num_atom_types,i,k,m,test,j,verbose; char **words; memory->sfree(params); @@ -388,14 +392,15 @@ void PairSMTBQ::read_file(char *file) verbose = 1; verbose = 0; - // open file on proc 0 + // open file on all processors FILE *fp; - fp = fopen( file, "r" ); + fp = force->open_potential(file); if ( fp == NULL ) { - fprintf( stderr, "error opening the force filed file! terminating...\n" ); + char str[128]; + sprintf(str,"Cannot open SMTBQ potential file %s",file); + error->one(FLERR,str); } - // read each line out of file, skipping blank lines or leading '#' // store line of params if all 3 element tags are in element list @@ -412,14 +417,14 @@ void PairSMTBQ::read_file(char *file) if (verbose) printf ("\n"); fgets(ptr,MAXLINE,fp); while (strchr(ptr,'#')) { - if (verbose) printf ("%s",ptr); - fgets(ptr,MAXLINE,fp); + if (verbose) printf ("%s",ptr); + fgets(ptr,MAXLINE,fp); } // Nombre d'atome different dans la structure // =============================================== - c = Tokenize( ptr, &words ); + Tokenize( ptr, &words ); num_atom_types = atoi(words[1]); if (verbose) printf (" %s %d\n", words[0], num_atom_types); @@ -430,7 +435,7 @@ void PairSMTBQ::read_file(char *file) 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; } + m = m + 1; } if (verbose) printf ("i %d, j %d, intype %d - nb pair %d\n",i,j,intype[i][j],m); } } @@ -440,10 +445,10 @@ void PairSMTBQ::read_file(char *file) if (nparams == maxparam) { maxparam += DELTA; params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), - "pair:params"); + "pair:params"); maxintparam += m; intparams = (Intparam *) memory->srealloc(intparams,(maxintparam+1)*sizeof(Intparam), - "pair:intparams"); + "pair:intparams"); } for (i=0; i < num_atom_types; i++) @@ -471,7 +476,7 @@ void PairSMTBQ::read_file(char *file) // Line 2 - Al fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + Tokenize( ptr, &words ); strcpy(params[i].nom , words[1]); params[i].sto = atof(words[2]); if (verbose) printf (" %s %s %f\n", words[0],params[i].nom,params[i].sto); @@ -479,7 +484,7 @@ void PairSMTBQ::read_file(char *file) //Line 3 - Charges fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + Tokenize( ptr, &words ); params[i].qform = atof(words[1]); params[i].masse = atof(words[2]); @@ -488,7 +493,7 @@ void PairSMTBQ::read_file(char *file) // Line 4 - Parametres QEq fgets( ptr, MAXLINE, fp); - c=Tokenize ( ptr, &words ); + Tokenize ( ptr, &words ); params[i].ne = atof(words[1]) ; params[i].chi = atof(words[2]) ; params[i].dj = atof(words[3]) ; @@ -496,7 +501,7 @@ void PairSMTBQ::read_file(char *file) if(strcmp(params[i].nom,"O")!=0){ params[i].R = atof(words[4]) ; if (verbose) printf(" %s %f %f %f %f\n",words[0],params[i].ne,params[i].chi, - params[i].dj,params[i].R); + params[i].dj,params[i].R); } else { if (verbose) printf(" %s %f %f %f\n",words[0],params[i].ne,params[i].chi,params[i].dj); } @@ -506,7 +511,7 @@ void PairSMTBQ::read_file(char *file) if(strcmp(params[i].nom,"O")==0){ fgets( ptr, MAXLINE, fp); - c=Tokenize ( ptr, &words ); + Tokenize ( ptr, &words ); coordOxBB= atof(words[1]) ; coordOxBulk= atof(words[2]) ; @@ -520,7 +525,7 @@ void PairSMTBQ::read_file(char *file) // Ligne 5 - Nombre d'etats partages fgets( ptr, MAXLINE, fp); - c=Tokenize ( ptr, &words ); + Tokenize ( ptr, &words ); params[i].n0 = atof(words[1]); if (verbose) printf(" %s %f\n",words[0],params[i].n0); @@ -534,7 +539,7 @@ void PairSMTBQ::read_file(char *file) reading the interaction's parameters ===================================================================== */ - m = 0; maxintsm = 0; // + m = 0; maxintsm = 0; // for (k=0 ; k<=maxintparam ; k++){intparams[k].intsm = 0;} // --------------------------------- for (k = 0; k < maxintparam; k++) { @@ -546,20 +551,20 @@ void PairSMTBQ::read_file(char *file) // Lecture des protagonistes fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + Tokenize( ptr, &words ); test = 0; 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); @@ -578,7 +583,7 @@ void PairSMTBQ::read_file(char *file) 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");} + error->all(FLERR,"the potential other than second_moment or buckingham have not treated here\n");} // On detemrine le type d'interaction @@ -589,29 +594,29 @@ void PairSMTBQ::read_file(char *file) 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"); } + 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 (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); + intparams[m].typepot,intparams[m].mode); fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + 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); + 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 ); + Tokenize( ptr, &words ); intparams[m].dc1 = atof(words[1]) ; intparams[m].dc2 = atof(words[2]) ; @@ -619,21 +624,21 @@ void PairSMTBQ::read_file(char *file) 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); - } + 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); + intparams[m].typepot); fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + 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); @@ -643,22 +648,22 @@ void PairSMTBQ::read_file(char *file) 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); + intparams[m].typepot); fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + 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 ); + 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); + intparams[m].bOO,intparams[m].r1OO,intparams[m].r2OO); } if (verbose) printf (" intsm %d \n",intparams[m].intsm); @@ -672,12 +677,12 @@ void PairSMTBQ::read_file(char *file) // Ligne 9 - rayon de coupure Electrostatique if (test == 0) { - fgets(ptr,MAXLINE,fp); - if (verbose) printf ("%s\n",ptr); + fgets(ptr,MAXLINE,fp); + if (verbose) printf ("%s\n",ptr); - fgets( ptr, MAXLINE, fp); + fgets( ptr, MAXLINE, fp); } - c= Tokenize( ptr, &words ); + Tokenize( ptr, &words ); for (i=0 ; i(kmax) ; if (verbose) printf (" kmax %d et ds %f\n",kmax,ds); /* ======================================================== */ @@ -700,12 +705,12 @@ void PairSMTBQ::read_file(char *file) if (verbose) printf ("%s",ptr); fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + Tokenize( ptr, &words ); Qstep = atoi(words[1]); - if (verbose) printf (" %s %d\n",words[0],Qstep); + if (verbose) printf (" %s " BIGINT_FORMAT "\n",words[0],Qstep); fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + Tokenize( ptr, &words ); loopmax = atoi(words[1]); precision = atof(words[2]); if (verbose) printf (" %s %d %f\n",words[0],loopmax,precision); @@ -716,22 +721,20 @@ void PairSMTBQ::read_file(char *file) if (verbose) printf ("%s",ptr); fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + Tokenize( ptr, &words ); r1Coord = atof(words[1]); r2Coord = atof(words[2]); if (verbose) printf (" %s %f %f\n",words[0],r1Coord,r2Coord); - - /* Mode for QInit============================================= */ fgets( ptr, MAXLINE, fp); if (verbose) printf ("%s",ptr); fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + Tokenize( ptr, &words ); strcpy( QInitMode , words[1] ); - if (strcmp(QInitMode,"true") == 0) { QOxInit= atof(words[2]); } - else { QOxInit = 0.0; } + if (strcmp(QInitMode,"true") == 0) QOxInit= atof(words[2]); + else QOxInit = 0.0; if (verbose) printf (" %s %s %f\n",words[0],QInitMode,QOxInit); @@ -741,32 +744,33 @@ void PairSMTBQ::read_file(char *file) if (verbose) printf ("%s",ptr); fgets( ptr, MAXLINE, fp); - c= Tokenize( ptr, &words ); + Tokenize( ptr, &words ); strcpy( QEqMode , words[1] ); if (verbose) printf (" %s %s\n",words[0],QEqMode); fgets( ptr, MAXLINE, fp); - if (strcmp(QEqMode,"BulkFromSlab") == 0) - { + if (strcmp(QEqMode,"BulkFromSlab") == 0) { + Tokenize( ptr, &words ); + zlim1QEq = atof(words[1]); + zlim2QEq = atof(words[2]); + if (verbose) printf (" %s %f %f\n",words[0],zlim1QEq,zlim2QEq); - c= Tokenize( ptr, &words ); - zlim1QEq = atof(words[1]); - zlim2QEq = atof(words[2]); - if (verbose) printf (" %s %f %f\n",words[0],zlim1QEq,zlim2QEq); + } else if (strcmp(QEqMode,"Surface") == 0) { + Tokenize( ptr, &words ); + zlim1QEq = atof(words[1]); + if (verbose) printf (" %s %f \n",words[0],zlim1QEq); - } - else if (strcmp(QEqMode,"Surface") == 0) - { - c= Tokenize( ptr, &words ); - zlim1QEq = atof(words[1]); - if (verbose) printf (" %s %f \n",words[0],zlim1QEq); - - } - else if (strcmp(QEqMode,"QEqAll") != 0 && - strcmp(QEqMode,"QEqAllParallel") != 0 && - strcmp(QEqMode,"Surface") != 0 ) - { error->all(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 zlim1all(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); - } + intparams[m].neig_cut = 1.2*intparams[m].r0; + if (strcmp(intparams[m].typepot,"second_moment") == 0 ) + if (verbose) printf (" Rc 1er voisin, typepot %s -> %f Ang\n", + intparams[m].typepot,intparams[m].neig_cut); + } } //A adapter au STO @@ -828,19 +830,17 @@ void PairSMTBQ::read_file(char *file) 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; + tagint itag,jtag; + int itype,jtype; int *ilist,*jlist,*numneigh,**firstneigh; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; @@ -849,10 +849,6 @@ void PairSMTBQ::compute(int eflag, int vflag) double **tmp,**tmpAll,*nmol; double dq,dqcov; - - int bavard; - - if (atom->nmax > nmax) { memory->destroy(ecov); memory->destroy(potmad); @@ -895,28 +891,32 @@ void PairSMTBQ::compute(int eflag, int vflag) double **x = atom->x; double **f = atom->f; double *q = atom->q; - int *tag = atom->tag; + +#ifdef LAMMPS_BIGBIG + long int *tag = atom->tag; +#else + tagint *tag = atom->tag; +#endif + 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); + const bigint step = update->ntimestep; + if (step == 0 || ((Qstep > 0) && (step % Qstep == 0))) Charge(); + // :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + // this is necessary to get sbcov or sbmet table in order to caclulate the covalent or metal bonding + if (Qstep == 0 || (step % Qstep != 0)) QForce_charge(0); // Charges Communication // ---------------------- - forward(q) ; reverse(q); + forward(q) ; // reverse(q); memory->create(nmol,nteam+1,"pair:nmol"); memory->create(tmp,nteam+1,7,"pair:tmp"); @@ -924,7 +924,7 @@ void PairSMTBQ::compute(int eflag, int vflag) for (i=0; i(nQEqall[i]); for (j=0; j<7; j++) { tmp[i][j] = 0.0; tmpAll[i][j] = 0.0; } } @@ -940,14 +940,10 @@ void PairSMTBQ::compute(int eflag, int vflag) 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]]; @@ -966,8 +962,8 @@ void PairSMTBQ::compute(int eflag, int vflag) tmp[gp][6] += iq*iq; -// self energy, only on i atom -// --------------------------- + // self energy, only on i atom + // --------------------------- Eself = self(¶ms[itype],iq); tmp[gp][0] += Eself; tmp[gp][2] += Eself; @@ -975,8 +971,8 @@ void PairSMTBQ::compute(int eflag, int vflag) if (evflag) ev_tally_full (i,0.0,2.0*Eself,0.0,0.0,0.0,0.0); -// N-body energy of i -// --------------------- + // N-body energy of i + // --------------------- dq = fabs(params[itype].qform) - fabs(iq); dqcov = dq*(2.0*ncov/params[itype].sto - dq); @@ -987,8 +983,8 @@ void PairSMTBQ::compute(int eflag, int vflag) if (evflag) ev_tally_full(i,0.0,2.0*ecov[i],0.0,0.0,0.0,0.0); -// Coulombian Interaction -// ----------------------- + // Coulombian Interaction + // ----------------------- evdwl = 2.0*Vself*iq*iq ; tmp[gp][1] += Vself*iq*iq; tmp[gp][2] += Vself*iq*iq; @@ -1001,13 +997,13 @@ void PairSMTBQ::compute(int eflag, int vflag) 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) { @@ -1017,7 +1013,7 @@ void PairSMTBQ::compute(int eflag, int vflag) 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 @@ -1031,9 +1027,9 @@ void PairSMTBQ::compute(int eflag, int vflag) rsq = delx*delx + dely*dely + delz*delz; -// --------------------------------- + // --------------------------------- if (sqrt(rsq) > cutmax) continue; -// --------------------------------- + // --------------------------------- // Coulombian Energy @@ -1056,106 +1052,106 @@ void PairSMTBQ::compute(int eflag, int vflag) if (evflag) - ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz); + 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; + 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; + // 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); + 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; + 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; + // 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); + 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; + // 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; + 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); + if (evflag) + ev_tally(i,j,nlocal,newton_pair,2.0*evdwl,0.0,fpair,delx,dely,delz); - evdwl = 0.0 ; fpair = 0.0; -// ----- ----- ----- ----- ----- ----- + evdwl = 0.0 ; fpair = 0.0; + // ----- ----- ----- ----- ----- ----- -// Attraction : force -// ------------------ - fpair = 0.0; - f_att(&intparams[m], i, j, rsq, fpair) ; + // 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; + 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); @@ -1173,31 +1169,31 @@ void PairSMTBQ::compute(int eflag, int vflag) 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) { + // :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + if (me == 0 && fmod(static_cast(step), Nevery) == 0.0 && strcmp(writeenerg,"true") == 0) { - ofstream fichierE; + 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 (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); + newton_pair,evflag,force->pair->tail_flag,vflag_fdotr); printf ("neighbor->includegroup %d\n",neighbor->includegroup); @@ -1207,8 +1203,8 @@ void PairSMTBQ::compute(int eflag, int vflag) printf ("Energie de la team %d -- %d atome -------(nteam = %d) \n ",gp,int(nmol[gp]),nteam); if (nmol[gp] == 0) { - printf (" ============================================ \n \n"); - continue; + printf (" ============================================ \n \n"); + continue; } printf ("Vself %f, Som q2 : %f, nmol %f\n",Vself,tmpAll[gp][6],nmol[gp]); // printf ("nmol %f\n",nmol[gp]); @@ -1234,7 +1230,6 @@ void PairSMTBQ::compute(int eflag, int vflag) memory->destroy(nmol); memory->destroy(tmp); memory->destroy(tmpAll); - } /* ---------------------------------------------------------------------- @@ -1373,7 +1368,7 @@ void PairSMTBQ::tabqeq() for (k=0; k < kmax+5; k++) //------------------------- { - s = double(k)*ds ; r = sqrt(s); + s = static_cast(k)*ds ; r = sqrt(s); if (k==0) r=10e-30; potqn[k] = 14.4*(erfc(alf*r)/r - mu) ; @@ -1432,107 +1427,107 @@ void PairSMTBQ::tabqeq() 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 ; + // -------------------------- + { + 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; + s = static_cast(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) ; + gammas(na,nb,za,zb,r,gam,dgam,dza,dzb, + d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ; - // --- Jij + // --- Jij - dij = 14.4 * (1.0/r - double(gam)); - ddij = 14.4 * (-1.0/(r*r) - double(dgam)) ; + dij = 14.4 * (1.0/r - static_cast(gam)); + ddij = 14.4 * (-1.0/(r*r) - static_cast(dgam)) ; - // Cutting Fonction + // Cutting Fonction - 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 (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) /square(nang); + } + if (r > rc) {dij = aCoeff *square(r- rc-nang) *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 (r > (rc+nang)) {dij = 0.0 ; ddij = 0.0;} - fafb[k][m] = potqn[k] - dij ; - if (k == 1) fafb[0][m] = fafb[k][m] ; + fafb[k][m] = potqn[k] - dij ; + if (k == 1) fafb[0][m] = fafb[k][m] ; - dfafb[k][m] = dpotqn[k] - ddij/r ; - } + 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); } + 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 ; + 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; + s = static_cast(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) ; + gammas(na,nb,za,zb,r,gam,dgam,dza,dzb, + d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ; - // --- Jij + // --- Jij - dij = 14.4 * (1.0/r - double(gam)); - ddij = 14.4 * (-1.0/(r*r) - double(dgam)) ; + dij = 14.4 * (1.0/r - static_cast(gam)); + ddij = 14.4 * (-1.0/(r*r) - static_cast(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 (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) /square(nang); + } + if (r > rc) {dij = aCoeff *square(r- rc-nang) *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 (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] ; + 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] ; + 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 ;} + dfafbTiOxSurf[k] = dpotqn[k] - ddij/r ;} - } + } } //for k @@ -1541,67 +1536,63 @@ void PairSMTBQ::tabqeq() // 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); } + 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 ; + 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; + s = static_cast(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) ; + gammas(na,nb,za,zb,r,gam,dgam,dza,dzb, + d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ; - // --- Jij + // --- Jij - dij = 14.4 * (1.0/r - double(gam)); - ddij = 14.4 * (-1.0/(r*r) - double(dgam)) ; + dij = 14.4 * (1.0/r - static_cast(gam)); + ddij = 14.4 * (-1.0/(r*r) - static_cast(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 (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) /square(nang); + } + if (r > rc) { + dij = aCoeff *square(r- rc-nang) *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 ; - } - } - + 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 @@ -1636,7 +1627,7 @@ void PairSMTBQ::potqeq(int i, int j, double qi, double qj, double rsq, r = rsq; sds = r/ds ; l = int(sds) ; - xi = sds - double(l) ; + xi = sds - static_cast(l) ; iCoord=coord[i]; @@ -1665,7 +1656,7 @@ void PairSMTBQ::potqeq(int i, int j, double qi, double qj, double rsq, 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 + // 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); @@ -1676,8 +1667,8 @@ void PairSMTBQ::potqeq(int i, int j, double qi, double qj, double rsq, 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))) ; + + (iIntfCoup2+jIntfCoup2)*((engBulk-engSurf)/(2*(coordOxBulk-coordOxSurf)) + - (engBB-engBulk)/(2*(coordOxBB-coordOxBulk))) ; // ---- Interpolation des forces @@ -1692,11 +1683,10 @@ void PairSMTBQ::potqeq(int i, int j, double qi, double qj, double rsq, 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))) ; + + (iIntfCoup2+jIntfCoup2)*((fforceBulk-fforceSurf)/(2*(coordOxBulk-coordOxSurf)) + - (fforceBB-fforceBulk)/(2*(coordOxBB-coordOxBulk))) ; - } - else{ // between metal and oxygen + } 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); @@ -1706,11 +1696,11 @@ void PairSMTBQ::potqeq(int i, int j, double qi, double qj, double rsq, 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; } + { iIntfCoup2=jIntfCoup2; + iCoord=jCoord; } eng = engBulk + (engBulk-engSurf)/(coordOxBulk-coordOxSurf) * iIntfCoup2 - + (engBB-engBulk)/(coordOxBB-coordOxBulk) * (iCoord-coordOxBulk-iIntfCoup2); + + (engBB-engBulk)/(coordOxBB-coordOxBulk) * (iCoord-coordOxBulk-iIntfCoup2); // ---- Forces Interpolation @@ -1730,18 +1720,13 @@ void PairSMTBQ::potqeq(int i, int j, double qi, double qj, double rsq, 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)); + - (engBulk-engSurf) *dIntfcoup2loc) + + 1/(coordOxBB-coordOxBulk) * ((fforceBB-fforceBulk)*(iCoord-coordOxBulk- iIntfCoup2) + - (engBB-engBulk) *(dcoordloc-dIntfcoup2loc)); } - - - } - - } /* -------------------------------------------------------------------- */ @@ -1758,8 +1743,8 @@ void PairSMTBQ::pot_ES (int i, int j, double rsq, double &eng) ==================================================================== */ int itype,jtype,l,m; - double r,t1,t2,sds,xi,engBulk,engSurf,dcoordloc,dcoupureloc; - double engBB, dIntfcoup2loc,iCoord,jCoord,iIntfCoup2,jIntfCoup2; + double r,t1,t2,sds,xi,engBulk,engSurf; + double engBB,iCoord,jCoord,iIntfCoup2,jIntfCoup2; int *type = atom->type; // int n = atom->ntypes; @@ -1770,7 +1755,7 @@ void PairSMTBQ::pot_ES (int i, int j, double rsq, double &eng) r = rsq; sds = r/ds ; l = int(sds) ; - xi = sds - double(l) ; + xi = sds - static_cast(l) ; iCoord=coord[i]; @@ -1801,8 +1786,8 @@ void PairSMTBQ::pot_ES (int i, int j, double rsq, double &eng) 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))) ; + + (iIntfCoup2+jIntfCoup2)*((engBulk-engSurf)/(2*(coordOxBulk-coordOxSurf)) + - (engBB-engBulk)/(2*(coordOxBB-coordOxBulk))) ; } else { // between metal and oxygen @@ -1815,12 +1800,12 @@ void PairSMTBQ::pot_ES (int i, int j, double rsq, double &eng) engBB = (t1 + (t2 - t1)*xi/2.0); if (jtype==0) { //the atom j is an oxygen - iIntfCoup2=jIntfCoup2; - iCoord=jCoord; - } + iIntfCoup2=jIntfCoup2; + iCoord=jCoord; + } eng = engBulk + (engBulk-engSurf)/(coordOxBulk-coordOxSurf)*iIntfCoup2 - + (engBB-engBulk)/(coordOxBB-coordOxBulk) * (iCoord-coordOxBulk-iIntfCoup2); + + (engBB-engBulk)/(coordOxBB-coordOxBulk) * (iCoord-coordOxBulk-iIntfCoup2); } @@ -1850,7 +1835,7 @@ void PairSMTBQ::pot_ES2 (int i, int j, double rsq, double &pot) r = rsq ; sds = r/ds ; l = int(sds) ; - xi = sds - double(l) ; + xi = sds - static_cast(l) ; // ---- Energies Interpolation @@ -1862,11 +1847,11 @@ void PairSMTBQ::pot_ES2 (int i, int j, double rsq, double &pot) } /* -------------------------------------------------------------------- - Oxygen-Oxygen Interaction + Oxygen-Oxygen Interaction -------------------------------------------------------------------- */ void PairSMTBQ::rep_OO(Intparam *intparam, double rsq, double &fforce, - int eflag, double &eng) + int eflag, double &eng) { double r,tmp_exp,tmp; double A = intparam->abuck ; @@ -1884,9 +1869,9 @@ void PairSMTBQ::rep_OO(Intparam *intparam, double rsq, double &fforce, void PairSMTBQ::Attr_OO(Intparam *intparam, double rsq, double &fforce, - int eflag, double &eng) + int eflag, double &eng) { - double r,tmp_exp,tmp; + double r,tmp_exp; double aOO = intparam->aOO ; double bOO = intparam->bOO ; double r1OO = intparam->r1OO ; @@ -1902,7 +1887,7 @@ void PairSMTBQ::Attr_OO(Intparam *intparam, double rsq, double &fforce, /* ---------------------------------------------------------------------- - covalente Interaction + covalente Interaction ----------------------------------------------------------------------*/ @@ -1933,63 +1918,63 @@ void PairSMTBQ::tabsm() 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)); + s = static_cast(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) - { + if (r <= rc1) + { - // -- Energy - tabsmb[k][sm] = Ksi*Ksi * tmpb ; - tabsmr[k][sm] = A * tmpr ; + // -- 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. */ + // -- 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 ; + dtabsmb[k][sm] = - 2.0 *Ksi*Ksi* q/rzero * tmpb /r; + dtabsmr[k][sm] = - A * p/rzero * tmpr/r ; - } // if + } // if - else if (r > rc1 && r <= rc2) - { + 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 ; + // -- 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. */ + // -- 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 ; + 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 - { + else + { - // -- Energie - tabsmb[k][sm] = 0.0; - tabsmr[k][sm] = 0.0; + // -- Energie + tabsmb[k][sm] = 0.0; + tabsmr[k][sm] = 0.0; - // -- Force - dtabsmb[k][sm] = 0.0; - dtabsmr[k][sm] = 0.0; + // -- Force + dtabsmb[k][sm] = 0.0; + dtabsmr[k][sm] = 0.0; - } + } @@ -2007,7 +1992,7 @@ void PairSMTBQ::tabsm() /* -------------------------------------------------------------- */ void PairSMTBQ::repulsive(Intparam *intparam, double rsq, int i, int j, - double &fforce, int eflag, double &eng) + double &fforce, int eflag, double &eng) { /* ================================================ @@ -2018,7 +2003,7 @@ void PairSMTBQ::repulsive(Intparam *intparam, double rsq, int i, int j, =================================================*/ int l; - double r,sds,xi,t1,t2,dt1,dt2,sweet,iq,jq; + double r,sds,xi,t1,t2,dt1,dt2,sweet; double rrcs = intparam->dc2; int sm = intparam->intsm; @@ -2030,7 +2015,7 @@ void PairSMTBQ::repulsive(Intparam *intparam, double rsq, int i, int j, if (sqrt(r) > rrcs) return ; sds = r/ds ; l = int(sds) ; - xi = sds - double(l) ; + xi = sds - static_cast(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) ; @@ -2057,11 +2042,11 @@ void PairSMTBQ::repulsive(Intparam *intparam, double rsq, int i, int j, void PairSMTBQ::attractive(Intparam *intparam, double rsq, - int eflag, int i, double iq, int j, double jq) + 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 sweet,mu; double rrcs = intparam->dc2; int *type = atom->type; @@ -2074,7 +2059,7 @@ void PairSMTBQ::attractive(Intparam *intparam, double rsq, sds = r/ds ; l = int(sds) ; - xi = sds - double(l) ; + xi = sds - static_cast(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) ; @@ -2082,24 +2067,24 @@ void PairSMTBQ::attractive(Intparam *intparam, double rsq, if (strcmp(intparam->mode,"oxide") == 0) { - mu = 0.5*(sqrt(params[1].sto) + sqrt(params[0].sto)); + 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); + // 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; + 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 (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 ; - } + sweet = 1.0; + sbmet[i] += (t1 + (t2 - t1)*xi/2.0) * sweet ; + } } @@ -2123,7 +2108,7 @@ void PairSMTBQ::f_att(Intparam *intparam, int i, int j,double rsq, double &fforc r = rsq; sds = r/ds ; l = int(sds) ; - xi = sds - double(l) ; + xi = sds - static_cast(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) ; @@ -2135,22 +2120,22 @@ void PairSMTBQ::f_att(Intparam *intparam, int i, int j,double rsq, double &fforc 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]) ) ; + + 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]) ) ; + + 1.0/sqrt(sbcov[j]*dqcovj + sbmet[j]) ) ; } } @@ -2172,9 +2157,9 @@ void PairSMTBQ::pot_COV(Param *param, int i, double &qforce) DQ = dq*(2.0*ncov/sto - dq); if (fabs(iq) < 1.0e-7 || fabs(sbcov[i]) < 1.0e-7) { - qforce = 0.0; } + qforce = 0.0; } else { - qforce = sign*sbcov[i]/sqrt(sbcov[i]*DQ + sbmet[i])*(ncov/sto - dq) ; + qforce = sign*sbcov[i]/sqrt(sbcov[i]*DQ + sbmet[i])*(ncov/sto - dq) ; } } @@ -2182,7 +2167,7 @@ void PairSMTBQ::pot_COV(Param *param, int i, double &qforce) /* ---------------------------------------------------------------------- */ double PairSMTBQ::potmet(Intparam *intparam, double rsq, - int i, double iq, int j, double jq) + int i, double iq, int j, double jq) { int l,itype,jtype; int *type = atom->type; @@ -2194,7 +2179,7 @@ double PairSMTBQ::potmet(Intparam *intparam, double rsq, r = rsq; sds = r/ds ; l = int(sds) ; - xi = sds - double(l) ; + xi = sds - static_cast(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) ; @@ -2207,14 +2192,14 @@ double PairSMTBQ::potmet(Intparam *intparam, double rsq, 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])) ); + + 1.0/(2.0*sqrt(sbcov[j]*dqcovj+sbmet[j])) ); return chi; } /* ---------------------------------------------------------------------- - Cutting Function + Cutting Function ----------------------------------------------------------------------*/ @@ -2242,7 +2227,7 @@ double PairSMTBQ::fcoupure(double r, double rep_dc1, double rep_dc2) } /* ---------------------------------------------------------------------- - Derivate of cutting function + Derivate of cutting function ----------------------------------------------------------------------*/ @@ -2301,7 +2286,7 @@ 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); + 6* c *(x*x - delta*delta)))/(16* delta*delta*delta); } @@ -2327,15 +2312,15 @@ double PairSMTBQ::Intfcoup2(double c, double x, double delta) /* --------------------------------------------------------------------- - Energy derivation respect charge Q + 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; + double xtmp,ytmp,ztmp; + double rsq; int *ilist,*jlist,*numneigh,**firstneigh; double iq,jq,fqi,fqj,fqij,fqij2,fqjj; int eflag; @@ -2351,89 +2336,85 @@ void PairSMTBQ::QForce_charge(int loop) 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) { - // ================== + // ================== + if (loop == 0) { + // ================== - for (ii = 0; ii < inum; ii ++) { - //-------------------------------- - i = ilist[ii]; - itype = map[type[i]]; + for (ii = 0; ii < inum; ii ++) { + //-------------------------------- + i = ilist[ii]; + itype = map[type[i]]; - gp = flag_QEq[i]; + gp = flag_QEq[i]; - sbcov[i] =coord[i]= sbmet[i] = 0.0; + 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]; + itype = map[type[i]]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + iq = q[i]; - // two-body interactions + // two-body interactions - jlist = firstneigh[i]; - jnum = numneigh[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (jj = 0; jj < jnum; jj++) { -// ------------------------------- - j = jlist[jj]; - j &= NEIGHMASK; + for (jj = 0; jj < jnum; jj++) { + // ------------------------------- + j = jlist[jj]; + j &= NEIGHMASK; - jtype = map[type[j]]; - jq = q[j]; + jtype = map[type[j]]; + jq = q[j]; - m = intype[itype][jtype]; - if (intparams[m].intsm == 0) continue ; + 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); + const double delx = x[j][0] - xtmp; + const double dely = x[j][1] - ytmp; + const double delz = x[j][2] - ztmp; + rsq = delx*delx + dely*dely + delz*delz; + + // Covalente charge forces - sbcov initialization + // ------------------------------------------------ + if (sqrt(rsq) > intparams[m].dc2) continue; + + attractive (&intparams[m],rsq,eflag,i,iq,j,jq); -// Covalente charge forces - sbcov initialization -// ------------------------------------------------ - if (sqrt(rsq) > intparams[m].dc2) continue; - - attractive (&intparams[m],rsq,eflag,i,iq,j,jq); + } // ---------------------------------------------- for jj - } // ---------------------------------------------- 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); - } // -------------------------------------------- ii + if (nteam == 0) return; //no oxide + if (Qstep == 0 || (step % Qstep != 0)) return; - // =============================================== - // 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) - // ======================= + // ======================= + } // End of If(loop) + // ======================= // =============================================== @@ -2455,12 +2436,12 @@ void PairSMTBQ::QForce_charge(int loop) // Madelung potential // -------------------- - potmad[i] += 2.0*Vself*iq ; + potmad[i] += 2.0*Vself*iq ; // charge force from self energy // ----------------------------- - fqi = qfo_self(¶ms[itype],iq); - potself[i] = fqi ; + fqi = qfo_self(¶ms[itype],iq); + potself[i] = fqi ; @@ -2473,29 +2454,29 @@ void PairSMTBQ::QForce_charge(int loop) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - jtype = map[type[j]]; - m = intype[itype][jtype]; - jq = q[j]; + jtype = map[type[j]]; + m = intype[itype][jtype]; + jq = q[j]; - delr[0] = x[j][0] - xtmp; - delr[1] = x[j][1] - ytmp; - delr[2] = x[j][2] - ztmp; - rsq = vec3_dot(delr,delr); + const double delx = x[j][0] - xtmp; + const double dely = x[j][1] - ytmp; + const double delz = x[j][2] - ztmp; + rsq = delx*delx + dely*dely + delz*delz; - // long range q-dependent - if (sqrt(rsq) > cutmax) continue; + // long range q-dependent + if (sqrt(rsq) > cutmax) continue; - // 1/r charge forces - // -------------------- - fqij = 0.0; -// pot_ES2 (i,j,rsq,fqij2); + // 1/r charge forces + // -------------------- + fqij = 0.0; + // pot_ES2 (i,j,rsq,fqij2); pot_ES (i,j,rsq,fqij); - potmad[i] += jq*fqij ; + potmad[i] += jq*fqij ; - } // ------ jj + } // ------ jj fqi = 0.0; pot_COV (¶ms[itype],i,fqi) ; @@ -2620,7 +2601,7 @@ void PairSMTBQ::Charge() i = ilist[ii]; itype = map[type[i]]; gp = flag_QEq[i]; -// if (gp == 0) continue; + // if (gp == 0) continue; if (itype == 0) q[i] = QOxInit ; if (itype > 0) q[i] = -QOxInit * nQEqaall[gp]/nQEqcall[gp]; @@ -2628,7 +2609,7 @@ void PairSMTBQ::Charge() } if (nteam == 0 || Qstep == 0) return; - if (fmod(double(step), double(Qstep)) != 0.0) return; + if (step % Qstep != 0) return; // -------------------------------------- // ---- @@ -2640,8 +2621,8 @@ void PairSMTBQ::Charge() 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 (" neutralite des charges %f\n qtotc %f qtota %f\n", + qtotll,qtotcll[gp]/nQEqcall[gp],qtotall[gp]/nQEqaall[gp]); printf (" ---------------------------- \n");} } @@ -2658,7 +2639,7 @@ void PairSMTBQ::Charge() // -------------------------------------------- for (iloop = 0; iloop < loopmax; iloop ++ ) { - // -------------------------------------------- + // -------------------------------------------- qtot = qtotll = Transf[3*cluster] = 0.0 ; for (gp=0; gp(nQEqall[i]); enegchk[i] = enegmax[i] = 0.0; } @@ -2729,9 +2710,9 @@ void PairSMTBQ::Charge() for (gp = 0; gp < nteam+1; gp++) { if(nQEqall[gp] !=0) { - enegchk[gp] = enegchkall[gp]/double(nQEqall[gp]); - enegmax[gp] = enegmaxall[gp]; - } + enegchk[gp] = enegchkall[gp]/static_cast(nQEqall[gp]); + enegmax[gp] = enegmaxall[gp]; + } } // ----------------------------------------------------- @@ -2744,7 +2725,7 @@ void PairSMTBQ::Charge() } for (gp = 1; gp < nteam+1; gp++) { m += end[gp] ; } - if (m == nteam) { break; } + if (m == nteam) { break; } // ----------------------------------------------------- // ----------------------------------------------------- @@ -2761,40 +2742,40 @@ void PairSMTBQ::Charge() // ======================================= // Charge Communication. // ======================================= - forward(q); reverse(q); + forward(q); // reverse(q); //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: // ========================================== // Ecriture des potentiels dans un fichier // ========================================== - if (strcmp(writepot,"true") == 0 && fmod(double(step), Neverypot) == 0.0) { + if (strcmp(writepot,"true") == 0 && fmod(static_cast(step), Neverypot) == 0.0) { - ofstream fichierpot("Electroneg_component.txt", ios::out | ios::trunc) ; + ofstream fichierpot("Electroneg_component.txt", ios::out | ios::trunc) ; - for (ii = 0; ii < inum; ii++) { + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; itype = map[type[i]]; gp = flag_QEq[i]; if (fichierpot) fichierpot<< setprecision(9) <(nQEqcall[i]) ; + TransfAll[i+2*cluster] /= static_cast(nQEqaall[i]) ;} sigmaa[i] = sigmac[i] = 0.0; } @@ -2816,8 +2797,8 @@ void PairSMTBQ::Charge() MPI_Allreduce(sigmac,sigmacll,nteam+1,MPI_DOUBLE,MPI_SUM,world); for (gp = 1; gp < nteam+1; gp++) { - sigmaall[gp] = sqrt(sigmaall[gp]/double(nQEqaall[gp])) ; - sigmacll[gp] = sqrt(sigmacll[gp]/double(nQEqcall[gp])) ; + sigmaall[gp] = sqrt(sigmaall[gp]/static_cast(nQEqaall[gp])) ; + sigmacll[gp] = sqrt(sigmacll[gp]/static_cast(nQEqcall[gp])) ; } @@ -2826,7 +2807,7 @@ void PairSMTBQ::Charge() for (gp = 0; gp < nteam+1; gp++) { printf (" -------------- Groupe %d -----------------\n",gp); printf (" qtotc %f(+- %f) qtota %f(+- %f)\n", - TransfAll[gp+cluster],sigmacll[gp],TransfAll[gp+2*cluster],sigmaall[gp]); + TransfAll[gp+cluster],sigmacll[gp],TransfAll[gp+2*cluster],sigmaall[gp]); printf (" Potentiel elec total : %f\n iloop %d, qtot %f\n",TransfAll[gp],iloop,TransfAll[3*cluster]); printf (" convergence : %f - %f\n",enegchk[gp],enegmax[gp]); } @@ -2858,9 +2839,9 @@ void PairSMTBQ::groupBulkFromSlab_QEq() i = ilist[ii]; ztmp = x[i][2]; if (ztmp>zlim1QEq && ztmp< zlim2QEq) - flag_QEq[i]=1; + flag_QEq[i]=1; else - flag_QEq[i]=0; + flag_QEq[i]=0; nteam=1; @@ -2901,9 +2882,9 @@ void PairSMTBQ::groupSurface_QEq() i = ilist[ii]; ztmp = x[i][2]; if (ztmp>zlim1QEq) - flag_QEq[i]=1; + flag_QEq[i]=1; else - flag_QEq[i]=0; + flag_QEq[i]=0; nteam=1; @@ -2986,22 +2967,22 @@ void PairSMTBQ::groupQEqAllParallel_QEq() jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++ ) - { - j = jlist[jj] ; jtype = map[type[j]]; - if (jtype == itype) continue; - m = intype[itype][jtype]; + { + j = jlist[jj] ; jtype = map[type[j]]; + if (jtype == itype) continue; + m = intype[itype][jtype]; - delr[0] = x[j][0] - xtmp; - delr[1] = x[j][1] - ytmp; - delr[2] = x[j][2] - ztmp; - rsq = vec3_dot(delr,delr); + 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) <= intparams[m].dc2) { - flag_QEq[i] = 1; flag_QEq[j] = 1; - } - } + if (sqrt(rsq) <= intparams[m].dc2) { + flag_QEq[i] = 1; flag_QEq[j] = 1; + } + } if (flag_QEq[i] == 1) { - QEq = 1; + QEq = 1; } } @@ -3063,84 +3044,84 @@ void PairSMTBQ::groupQEqAllParallel_QEq() // ---------------------- 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; + 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]; + 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]; + // 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 (jtype == 0 && flag_QEq[j] == 0) continue; - if (flag_gp[me][j] == ngp) 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); + 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) <= cutmax) { - flag_QEq[j] = 1; //Entre dans le schema QEq + flag_QEq[j] = 1; //Entre dans le schema QEq - // :::::::::::::::::::: Meeting of two group in the same proc ::::::::::::::::::::: + // :::::::::::::::::::: 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); + 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); + // 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; - } + 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; - } + 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; - } + 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; - } - // :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 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 + flag_gp[me][j] = ngp; + if (j < nlocal) + { + tab_gp[ngp][nelt[ngp]] = j; + nelt[ngp] += 1; + } + } + } // for j } // for k } // for ii @@ -3148,7 +3129,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq() // Groups communication // OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO for (i = 0; i < nproc; i++) { - forward_int(flag_gp[i]); reverse_int(flag_gp[i]); + forward_int(flag_gp[i]); //reverse_int(flag_gp[i]); } // --- @@ -3170,73 +3151,73 @@ void PairSMTBQ::groupQEqAllParallel_QEq() jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++ ) - { - j = jlist[jj] ; jtype = map[type[j]]; - if (jtype != 0) continue; + { + 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; + 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); + 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) { - // ---------------------------------------- + // ---------------------------------------- + if (sqrt(rsq) <= cutmax) { + // if (sqrt(rsq) <= intparams[m].dc2) { + // ---------------------------------------- - flag_QEq[i] = 1; igp = flag_gp[me][j]; + 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] == 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]); + 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]; } + 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 < 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; - } + 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; - } + } + 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; - } + 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 + } // voisin j } // atom i /* ================================================== - Group Communication between proc for unification + Group Communication between proc for unification ================================================== */ for (i = 0; i < nproc; i++) { - forward_int(flag_gp[i]); reverse_int(flag_gp[i]); + forward_int(flag_gp[i]);// reverse_int(flag_gp[i]); } // =============== End of COMM ================= @@ -3253,7 +3234,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq() if (tabtemp[m][k] != 0) continue; if (flag_gp[k][i] != 0) { - tabtemp[m][k] = 10*k + flag_gp[k][i]; + tabtemp[m][k] = 10*k + flag_gp[k][i]; } } @@ -3266,7 +3247,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq() nteam = 0; iproc = 0; for (igp = 0; igp < 10*nproc; igp++) { if (Allgptmp[igp] == 0) continue; - iproc = int(double(igp)/10.0); + iproc = int(static_cast(igp)/10.0); ngp = igp - 10*iproc; if (nteam == 0) { @@ -3277,17 +3258,17 @@ void PairSMTBQ::groupQEqAllParallel_QEq() } 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; - } } + 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; + nteam += 1; + team_elt[nteam][iproc] = 0; + team_QEq[nteam][iproc][team_elt[nteam][iproc]] = ngp; + team_elt[nteam][iproc] += 1; } } // ------- @@ -3297,31 +3278,31 @@ void PairSMTBQ::groupQEqAllParallel_QEq() 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 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; + // 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; + 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; - } + // 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 + 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 @@ -3339,8 +3320,8 @@ void PairSMTBQ::groupQEqAllParallel_QEq() 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; + flag_QEq[i] = j; m = 1; + break; } } if (m == 1) break; @@ -3404,7 +3385,7 @@ void PairSMTBQ::Init_charge(int *nQEq, int *nQEqa, int *nQEqc) 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]); + if (itype != 0) q[i] = 1.96 * static_cast(nQEqaall[gp]) / static_cast(nQEqcall[gp]); } tot += q[i]; } @@ -3424,9 +3405,9 @@ int PairSMTBQ::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, in 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]); + // if (j < 3) printf ("[%d] %d pfc %d %d buf_send = %f \n",me,n,i,m-1,buf[m-1]); } - return 1; + return m; } /* ---------------------------------------------------------------------- */ @@ -3439,7 +3420,7 @@ void PairSMTBQ::unpack_forward_comm(int n, int first, double *buf) last = first + n ; for (i = first; i < last; i++) { tab_comm[i] = buf[m++]; -// if (inlocal; int nghost = atom->nghost; - for (i=0; i(tab[i]);} comm->forward_comm_pair(this); @@ -3527,7 +3508,7 @@ void PairSMTBQ::reverse_int(int *tab) int nlocal = atom->nlocal; int nghost = atom->nghost; - for (i=0; i(tab[i]);} comm->reverse_comm_pair(this); @@ -3595,7 +3576,7 @@ void PairSMTBQ::CheckEnergyVSForce() double za,zb,gam,dgam,dza,dzb, d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2,na,nb; int *type = atom->type; - char *NameFile; + const char *NameFile; i=0; @@ -3625,25 +3606,25 @@ void PairSMTBQ::CheckEnergyVSForce() { if (iCoord==1) - {coord[i]=2.2; - coord[j]=2.1; - NameFile="energyandForceOxOxUnderCoord.txt"; - } + {coord[i]=2.2; + coord[j]=2.1; + NameFile=(const char *)"energyandForceOxOxUnderCoord.txt"; + } if (iCoord==2) - {coord[i]=coordOxBulk; - coord[j]=coordOxBulk; - NameFile="energyandForceOxOxCoordBulk.txt"; - } + {coord[i]=coordOxBulk; + coord[j]=coordOxBulk; + NameFile=(const char *)"energyandForceOxOxCoordBulk.txt"; + } if (iCoord==3) - {coord[i]=3.8; - coord[j]=4; - NameFile="energyandForceOxOxOverCoord.txt"; - } + {coord[i]=3.8; + coord[j]=4; + NameFile=(const char *)"energyandForceOxOxOverCoord.txt"; + } if (iCoord==4) - {coord[i]=2.2; - coord[j]=3.5; - NameFile="energyandForceOxOxUnderOverCoord.txt"; - } + {coord[i]=2.2; + coord[j]=3.5; + NameFile=(const char *)"energyandForceOxOxUnderOverCoord.txt"; + } ofstream fichierOxOx(NameFile, ios::out | ios::trunc) ; @@ -3651,56 +3632,56 @@ void PairSMTBQ::CheckEnergyVSForce() 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; + 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; + 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 ; + 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) ; gammas(na,nb,za,zb,r,gam,dgam,dza,dzb, - d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ; + d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ; - sds = rsq/ds ; l = int(sds) ; - xi = sds - double(l) ; + sds = rsq/ds ; l = int(sds) ; + xi = sds - static_cast(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 = 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 = 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 = 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 = 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 = 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 ; + 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) <(na)*drtrm*ss; + dza=rtrm*dzeta1; + d2zaa=rtrm*d2zeta11; + d2zra=rtrm*d2zeta1r+static_cast(na)*drtrm*dzeta1; + d2gamr2=d2gamr2+0.5*static_cast(na*(na2-1))*d2rtrm*ss + 2.0*static_cast(na)*drtrm*deriv+rtrm*deriv2; -// Sum over 2*nb + // Sum over 2*nb - rtrm=drtrm; - drtrm=d2rtrm; - ztrm=0.5/(zb*float(nb2)); + rtrm=drtrm; + drtrm=d2rtrm; + ztrm=0.5/(zb*static_cast(nb2)); - for (i = nb2; i >= 1; i--) { - rtrm=rtrm*halfr; - drtrm=drtrm*halfr; - ztrm=ztrm*2.0*zb; - ctrm=ztrm/factorial(nb2-i); + for (i = nb2; i >= 1; i--) { + rtrm=rtrm*halfr; + drtrm=drtrm*halfr; + ztrm=ztrm*2.0*zb; + ctrm=ztrm/factorial(static_cast(nb2-i)); - css(ss,na2-1,nb2-i,z2ra,z2rb,r,deriv,dzeta1,dzeta2, - d2zeta11,d2zeta12,d2zeta22,d2zeta1r,d2zeta2r,deriv2); + 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); - } + trm1=static_cast(i)*ctrm; + trm2=trm1*rtrm; + gam=gam-trm2*ss; + trm3=trm1*static_cast(na2+nb2-i)*drtrm; + dgam=dgam-trm2*deriv-0.5*trm3*ss; + d2gamr2=d2gamr2-trm2*deriv2-trm3*deriv-0.5*trm3*static_cast(na2+nb2-i-1)*ss/r; + dza=dza-trm2*dzeta1; + dzb=dzb-(trm2/zb)*((static_cast(nb2-i))*ss+zb*dzeta2); + d2zaa=d2zaa-trm2*d2zeta11; + d2zab=d2zab-(trm2/zb)*((static_cast(nb2-i))*dzeta1+zb*d2zeta12); + d2zbb=d2zbb-(trm2/zb)*(2.0*(static_cast(nb2-i))*dzeta2+zb*d2zeta22 + + (static_cast((nb2-i-1)*(nb2-i))*ss/zb)); + d2zra=d2zra-trm2*d2zeta1r-0.5*trm3*dzeta1; + d2zrb=d2zrb-(trm2/zb)*((static_cast(nb2-i))*deriv+zb*d2zeta2r) - + 0.5*(trm3/zb)*((static_cast(nb2-i))*ss+zb*dzeta2); + } -// Multiply by coefficients + // 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; + trm3=powint(2.0*za,na2+1)/factorial(static_cast(na2)); + gam=gam*trm3; + dgam=dgam*trm3; + rfct1=((static_cast(na2+1))/za); + rgam1=rfct1*gam; + dza=dza*trm3; + rdza1=2.0*rfct1*dza; + dza=dza+rgam1; + dzb=dzb*trm3; + rgam2=rgam1*static_cast(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 + 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. + // 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 + 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 - ------------------------------------------------------------------- */ + 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; + 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]; @@ -3956,163 +3937,163 @@ void PairSMTBQ::css(double &s, double nn1, double nn2, double alpha, double beta memory->create(b,31,"pair:a"); -// Set up factorials - stored as factorial(n) in location(n+1) + // 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; + 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 + // 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; + 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 + // 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; - } + 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 + // 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); + 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 + // 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; - } + 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 + // 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; + 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); + memory->destroy(a); + memory->destroy(b); - return; + return; } /* ------------------------------------------------------------------------------- - coeffs + coeffs ------------------------------------------------------------------------------- */ - double PairSMTBQ::coeffs(int na, int nb, int k) - { - // implicit real (a-h,o-z) - // common /fctrl/ fct(30) +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; + int il,je,ia,i,j,ie,l; + double coeffs; -// Statement function -// binm(n,i)=fct(n+1)/(fct(n-i+1)*fct(i+1)); + // 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; - } + 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)*powint(-1.,j); + } + return coeffs; +} // ============================================ @@ -4122,120 +4103,96 @@ double PairSMTBQ::binm(int n, int i) } /* --------------------------------------------------------------------------------- - Caintgs + Caintgs -------------------------------------------------------------------------------- */ - void PairSMTBQ::caintgs (double x, int k, double *a) - { -// implicit real (a-h,o-z) -// dimension a(30) - int i; - double cste,rx; +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; - } + cste=exp(-x); + rx=1.0/x; + a[1]=cste*rx; + for (i = 1; i <= k; i++) { + a[i+1]=(a[i]*static_cast(i)+cste)*rx; + } + return; +} /* ----------------------------------------------------------------------------------- - Cbintgs + 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) + // 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; + int i0,m,last,i; + double absx,expx,expmx,ytrm,y,rx; - i0=0; - absx=fabs(x); + 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; + 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; + 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]=(static_cast(i)*b[i]+ powint(-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 = powint(-x,m-1)*(1. - powint(-1.,m+i+1))/(fct[m+1]*static_cast(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.-powint(-1.,i+1))/static_cast(i+1); + } + g190: + return; } - /* ============================== This is the END... ================================== */ diff --git a/src/USER-SMTBQ/pair_smtbq.h b/src/USER-SMTBQ/pair_smtbq.h index 2a1dc3bb8a..49fbe8768b 100644 --- a/src/USER-SMTBQ/pair_smtbq.h +++ b/src/USER-SMTBQ/pair_smtbq.h @@ -1,197 +1,195 @@ -/* ---------------------------------------------------------------------- - 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 const T& min ( const T& a, const T& b ) { + return !(b