1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 /* ----------------------------------------------------------------------
16 The SMTBQ code has been developed with the financial support of CNRS and
17 of the Regional Council of Burgundy (Convention n¡ 2010-9201AAO037S03129)
18
19 Copyright (2015)
20 Universite de Bourgogne : Nicolas SALLES, Olivier POLITANO
21 Universite de Paris-Sud Orsay : R. Tetot
22 Aalto University (Finland) : E. Maras
23
24 Please cite the related publication:
25 N. Salles, O. Politano, E. Amzallag and R. Tetot,
26 Comput. Mater. Sci., 111 (2016) 181-189
27
28 Contact : lammps@u-bourgogne.fr
29
30 This program is free software; you can redistribute it and/or
31 modify it under the terms of the GNU General Public License as
32 published by the Free Software Foundation; either version 2 of
33 the License, or (at your option) any later version.
34
35 This program is distributed in the hope that it will be useful,
36 but WITHOUT ANY WARRANTY; without even the implied warranty of
37 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
38 See the GNU General Public License for more details:
39 <https://www.gnu.org/licenses/>.
40 ------------------------------------------------------------------------- */
41
42 #include "pair_smtbq.h"
43
44 #include "atom.h"
45 #include "comm.h"
46 #include "error.h"
47 #include "force.h"
48 #include "math_const.h"
49 #include "math_extra.h"
50 #include "math_special.h"
51 #include "memory.h"
52 #include "neigh_list.h"
53 #include "neigh_request.h"
54 #include "neighbor.h"
55 #include "update.h"
56
57 #include <cmath>
58 #include <cstring>
59
60 #include <algorithm>
61 #include <fstream>
62 #include <iomanip>
63
64 using namespace std;
65
66 using namespace LAMMPS_NS;
67 using namespace MathConst;
68 using namespace MathExtra;
69 using namespace MathSpecial;
70
71 #define MAXLINE 2048
72 #define MAXTOKENS 2048
73 #define DELTA 4
74 #define PGDELTA 1
75 #define MAXNEIGH 24
76
77 /* ---------------------------------------------------------------------- */
78
PairSMTBQ(LAMMPS * lmp)79 PairSMTBQ::PairSMTBQ(LAMMPS *lmp) : Pair(lmp)
80 {
81 MPI_Comm_rank(world,&me);
82 MPI_Comm_size(world,&nproc);
83
84 single_enable = 0;
85 restartinfo = 0;
86 one_coeff = 1;
87
88 nmax = 0;
89 rmin = 0.0;
90 dr = 0.0;
91 ds = 0.0;
92 kmax = 0;
93
94 params = nullptr;
95 intparams = nullptr;
96
97 intype = nullptr;
98 coultype = nullptr;
99 fafb = nullptr;
100 dfafb = nullptr;
101 potqn = nullptr;
102 dpotqn = nullptr;
103 Vself = 0.0;
104 tabsmb = nullptr;
105 tabsmr = nullptr;
106 dtabsmb = nullptr;
107 dtabsmr = nullptr;
108
109 sbcov = nullptr;
110 coord = nullptr;
111 sbmet = nullptr;
112 ecov = nullptr;
113
114 potmad = nullptr;
115 potself = nullptr;
116 potcov = nullptr;
117 qf = nullptr;
118 q1 = nullptr;
119 q2 = nullptr;
120 tab_comm = nullptr;
121
122 nvsm = nullptr;
123 vsm = nullptr;
124 flag_QEq = nullptr;
125 nQEqaall = nullptr;
126 nQEqcall = nullptr;
127 nQEqall = nullptr;
128 nteam = 0;
129 cluster = 0;
130
131 Nevery = 0.0;
132 Neverypot = 0.0;
133
134 fct = nullptr;
135
136 maxpage = 0;
137
138 // set comm size needed by this Pair
139
140 comm_forward = 1;
141 comm_reverse = 1;
142 }
143
144 /* ----------------------------------------------------------------------
145 check if allocated, since class can be destructed when incomplete
146 ------------------------------------------------------------------------- */
147
~PairSMTBQ()148 PairSMTBQ::~PairSMTBQ()
149 {
150 int i;
151 if (elements) {
152 for (i = 0; i < atom->ntypes ; i++ ) free( params[i].nom);
153 for (i = 1; i <= maxintparam ; i++ ) free( intparams[i].typepot);
154 for (i = 1; i <= maxintparam ; i++ ) free( intparams[i].mode);
155 }
156
157 free(QEqMode);
158 free(QInitMode);
159 free(writepot);
160 free(writeenerg);
161 free(Bavard);
162
163 memory->sfree(params);
164 memory->sfree(intparams);
165
166 memory->destroy(intype);
167 memory->destroy(coultype);
168 memory->destroy(fafb);
169 memory->destroy(dfafb);
170 memory->destroy(potqn);
171 memory->destroy(dpotqn);
172
173 memory->destroy(fafbOxOxSurf);
174 memory->destroy(dfafbOxOxSurf);
175 memory->destroy(fafbTiOxSurf);
176 memory->destroy(dfafbTiOxSurf);
177
178 memory->destroy(fafbOxOxBB);
179 memory->destroy(dfafbOxOxBB);
180 memory->destroy(fafbTiOxBB);
181 memory->destroy(dfafbTiOxBB);
182
183 memory->destroy(ecov);
184 memory->destroy(sbcov);
185 memory->destroy(coord);
186 memory->destroy(sbmet);
187 memory->destroy(tabsmb);
188 memory->destroy(tabsmr);
189 memory->destroy(dtabsmb);
190 memory->destroy(dtabsmr);
191
192 memory->destroy(potmad);
193 memory->destroy(potself);
194 memory->destroy(potcov);
195
196 memory->destroy(nvsm);
197 memory->destroy(vsm);;
198 memory->destroy(flag_QEq);
199
200 memory->destroy(nQEqall);
201 memory->destroy(nQEqcall);
202 memory->destroy(nQEqaall);
203
204 memory->destroy(qf);
205 memory->destroy(q1);
206 memory->destroy(q2);
207 memory->destroy(tab_comm);
208
209 if (allocated) {
210 memory->destroy(setflag);
211 memory->destroy(cutsq);
212 delete [] esm;
213 }
214
215 memory->destroy(fct);
216 }
217
218 /* ---------------------------------------------------------------------- */
219
allocate()220 void PairSMTBQ::allocate()
221 {
222 allocated = 1;
223 int n = atom->ntypes;
224
225 memory->create(setflag,n+1,n+1,"pair:setflag");
226 memory->create(cutsq,n+1,n+1,"pair:cutsq");
227
228 map = new int[n+1];
229 esm = new double[n];
230 }
231
232 /* ----------------------------------------------------------------------
233 global settings
234 ------------------------------------------------------------------------- */
235
settings(int narg,char **)236 void PairSMTBQ::settings(int narg, char ** /* arg */)
237 {
238 if (narg > 0) error->all(FLERR,"Illegal pair_style command");
239 }
240
241 /* ----------------------------------------------------------------------
242 set coeffs for one or more type pairs
243 ------------------------------------------------------------------------- */
244
coeff(int narg,char ** arg)245 void PairSMTBQ::coeff(int narg, char **arg)
246 {
247 if (!allocated) allocate();
248
249 if (utils::strmatch(force->pair_style,"^hybrid"))
250 error->all(FLERR,"Pair style SMTBQ is not compatible with hybrid styles");
251
252 map_element2type(narg-3,arg+3);
253
254 // read potential file and initialize potential parameters
255
256 read_file(arg[2]);
257
258 // generate Coulomb 1/r energy look-up table
259
260 if (comm->me == 0 && screen)
261 fputs("Pair SMTBQ: generating Coulomb integral lookup table ...\n",screen);
262
263 tabqeq();
264
265 // ------------
266
267 if (comm->me == 0 && screen)
268 fputs(" generating Second Moment integral lookup table ...\n",screen);
269
270 tabsm();
271 }
272
273 /* ----------------------------------------------------------------------
274 init specific to this pair style
275 ------------------------------------------------------------------------- */
276
init_style()277 void PairSMTBQ::init_style()
278 {
279 if (atom->tag_enable == 0)
280 error->all(FLERR,"Pair style SMTBQ requires atom IDs");
281 if (force->newton_pair == 0)
282 error->all(FLERR,"Pair style SMTBQ requires newton pair on");
283 if (!atom->q_flag)
284 error->all(FLERR,"Pair style SMTBQ requires atom attribute q");
285
286
287 // need a full neighbor list
288
289 int irequest = neighbor->request(this);
290 neighbor->requests[irequest]->half = 0;
291 neighbor->requests[irequest]->full = 1;
292
293 pgsize = neighbor->pgsize;
294 oneatom = neighbor->oneatom;
295 // if (maxpage == 0) add_pages();
296
297 }
298
299 /* ---------------------------------------------------------------------- */
300
init_one(int i,int j)301 double PairSMTBQ::init_one(int i, int j)
302 {
303 if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
304 return cutmax;
305 }
306
307 /* ----------------------------------------------------------------------
308 ---------------------------------------------------------------------- */
309
read_file(char * file)310 void PairSMTBQ::read_file(char *file)
311 {
312 int num_atom_types,i,k,m,test,j,verbose;
313 char **words;
314
315 memory->sfree(params);
316 params = nullptr;
317 memory->sfree(intparams);
318 intparams = nullptr;
319 nparams = 0;
320 maxparam = 0;
321 maxintparam = 0;
322
323 verbose = 1;
324 verbose = 0;
325
326 coordOxBB = 0.0;
327 coordOxBulk = 0.0;
328 coordOxSurf = 0.0;
329 ROxBB = 0.0;
330 ROxSurf = 0.0;
331
332 // open file on all processors
333 FILE *fp;
334 fp = utils::open_potential(file,lmp,nullptr);
335 if (fp == nullptr) {
336 char str[128];
337 snprintf(str,128,"Cannot open SMTBQ potential file %s",file);
338 error->one(FLERR,str);
339 }
340
341 // read each line out of file, skipping blank lines or leading '#'
342 // store line of params if all 3 element tags are in element list
343
344 char *ptr;
345
346 ptr = (char*) malloc(sizeof(char)*MAXLINE);
347 words = (char**) malloc(sizeof(char*)*MAXTOKENS);
348 for (i=0; i < MAXTOKENS; i++)
349 words[i] = (char*) malloc(sizeof(char)*MAXTOKENS);
350
351
352 /* strip comment, skip line if blank */
353
354 if (verbose) printf ("\n");
355 fgets(ptr,MAXLINE,fp);
356 while (strchr(ptr,'#')) {
357 if (verbose) printf ("%s",ptr);
358 fgets(ptr,MAXLINE,fp);
359 }
360
361
362 // Nombre d'atome different dans la structure
363 // ===============================================
364 Tokenize( ptr, &words );
365 num_atom_types = atoi(words[1]);
366 if (verbose) printf (" %s %d\n", words[0], num_atom_types);
367
368 memory->create(intype,num_atom_types,num_atom_types,"pair:intype");
369
370 m = 0;
371 for (i = 0; i < num_atom_types; i++) {
372 for (j = 0; j < num_atom_types; j++) {
373 if (j < i) { intype[i][j] = intype[j][i];}
374 else { intype[i][j] = 0;
375 m = m + 1; }
376 if (verbose) printf ("i %d, j %d, intype %d - nb pair %d\n",i,j,intype[i][j],m);
377 }
378 }
379
380 // load up parameter settings and error check their values
381
382 nparams = maxparam = num_atom_types;
383 params = (Param *) memory->create(params,maxparam*sizeof(Param),
384 "pair:params");
385 maxintparam = m;
386 intparams = (Intparam *) memory->create(intparams,(maxintparam+1)*sizeof(Intparam),
387 "pair:intparams");
388
389 for (i=0; i < num_atom_types; i++)
390 params[i].nom = (char*) malloc(sizeof(char)*3);
391
392 for (i=1; i <= maxintparam; i++)
393 intparams[i].typepot = (char*) malloc(sizeof(char)*15);
394
395 for (i=1; i <= maxintparam; i++)
396 intparams[i].mode = (char*) malloc(sizeof(char)*6);
397
398 QEqMode = (char*) malloc(sizeof(char)*19);
399 Bavard = (char*) malloc(sizeof(char)*6);
400 QInitMode = (char*) malloc(sizeof(char)*19);
401 writepot = (char*) malloc(sizeof(char)*6);
402 writeenerg = (char*) malloc(sizeof(char)*6);
403
404
405 // Little loop for ion's parameters
406 // ================================================
407 for (i=0; i<num_atom_types; i++) {
408
409 fgets(ptr,MAXLINE,fp); if (verbose) printf ("%s",ptr);
410
411 // Line 2 - Al
412
413 fgets( ptr, MAXLINE, fp);
414 Tokenize( ptr, &words );
415 strcpy(params[i].nom , words[1]);
416 params[i].sto = atof(words[2]);
417 if (verbose) printf (" %s %s %f\n", words[0],params[i].nom,params[i].sto);
418
419 //Line 3 - Charges
420
421 fgets( ptr, MAXLINE, fp);
422 Tokenize( ptr, &words );
423
424 params[i].qform = atof(words[1]);
425 params[i].masse = atof(words[2]);
426 if (verbose) printf (" %s %f %f \n", words[0],params[i].qform, params[i].masse);
427
428 // Line 4 - Parametres QEq
429
430 fgets( ptr, MAXLINE, fp);
431 Tokenize ( ptr, &words );
432 params[i].ne = atof(words[1]) ;
433 params[i].chi = atof(words[2]) ;
434 params[i].dj = atof(words[3]) ;
435
436 if (strcmp(params[i].nom,"O")!=0) {
437 params[i].R = atof(words[4]) ;
438 if (verbose) printf(" %s %f %f %f %f\n",words[0],params[i].ne,params[i].chi,
439 params[i].dj,params[i].R);
440 } else {
441 if (verbose) printf(" %s %f %f %f\n",words[0],params[i].ne,params[i].chi,params[i].dj);
442 }
443
444
445 // Line 4bis - Coordinance et rayon pour Ox
446 if (strcmp(params[i].nom,"O")==0) {
447
448 fgets( ptr, MAXLINE, fp);
449 Tokenize ( ptr, &words );
450
451 coordOxBB= atof(words[1]) ;
452 coordOxBulk= atof(words[2]) ;
453 coordOxSurf= atof(words[3]) ;
454 ROxBB = atof(words[4]) ;
455 params[i].R = atof(words[5]) ;
456 ROxSurf = atof(words[6]) ;
457 if (verbose) printf(" %s %f %f %f %f %f %f\n",words[0],coordOxBB,coordOxBulk,coordOxSurf,ROxBB,params[i].R,ROxSurf);
458 }
459
460 // Ligne 5 - Nombre d'etats partages
461
462 fgets( ptr, MAXLINE, fp);
463 Tokenize ( ptr, &words );
464 params[i].n0 = atof(words[1]);
465 if (verbose) printf(" %s %f\n",words[0],params[i].n0);
466
467 // Parametres de Slater
468 params[i].dzeta = (2.0*params[i].ne + 1.0)/(4.0*params[i].R);
469 if (verbose) printf (" Parametre dzeta (Slater) : %f\n",params[i].dzeta);
470
471 } // Fin elements i
472
473 /* =====================================================================
474 reading the interaction's parameters
475 ===================================================================== */
476
477 m = 0; maxintsm = 0; //
478 for (k=0 ; k<=maxintparam ; k++) {intparams[k].intsm = 0;}
479 // ---------------------------------
480 for (k = 0; k < maxintparam; k++) {
481 // ---------------------------------
482 m += 1;
483
484 // Ligne 5 - parametre des potentials
485 fgets(ptr,MAXLINE,fp); if (verbose) printf ("%s",ptr);
486
487 // Lecture des protagonistes
488 fgets( ptr, MAXLINE, fp);
489 Tokenize( ptr, &words );
490
491 test = 0;
492 for (i = 0; i <num_atom_types; i++)
493 {
494 if (strcmp(params[i].nom,words[1])==0) break;
495 if (i == num_atom_types - 1) test = 1;
496 }
497 // if (test == 0) printf (" on a %s -> %d = %s\n",words[1],i,params[i].nom);
498
499 for (j = 0; j <num_atom_types; j++)
500 {
501 if (strcmp(params[j].nom,words[2])==0) break;
502 if (j == num_atom_types - 1) test = 1;
503 }
504 // if (test == 0) printf (" on a %s -> %d = %s\n",words[2],j,params[j].nom);
505
506
507 if (test == 1) {
508 if (verbose) printf ("========== fin des interaction ==========\n");
509 break ; }
510
511
512 intype[i][j] = m;
513 intype[j][i] = intype[i][j];
514 strcpy( intparams[m].typepot , words[3] );
515 intparams[m].intsm = 0;
516 if (verbose) printf (" itype %d jtype %d - intype %d\n",i,j,intype[i][j]);
517
518 if (strcmp(intparams[m].typepot,"second_moment") !=0 &&
519 strcmp(intparams[m].typepot,"buck") != 0 &&
520 strcmp(intparams[m].typepot,"buckPlusAttr") != 0) {
521 error->all(FLERR,"the potential other than second_moment or buckingham have not treated here\n");}
522
523
524 // On detemrine le type d'interaction
525 // -----------------------------------
526 if (strcmp(intparams[m].typepot,"second_moment") == 0) {
527 maxintsm += 1;
528 strcpy( intparams[m].mode , words[4] );
529 intparams[m].intsm = maxintsm;
530
531 if (strcmp(intparams[m].mode,"oxide") != 0 &&
532 strcmp(intparams[m].mode,"metal") != 0) {
533 error->all(FLERR,"needs mode to second moment interaction : oxide or metal"); }
534
535 // if (strcmp(intparams[m].mode,"oxide") == 0)
536 // intparams[m].ncov = min((params[i].sto)*(params[i].n0),(params[j].sto)*(params[j].n0));
537
538 if (verbose) printf(" %s %s %s %s %s \n",words[0],words[1],words[2],
539 intparams[m].typepot,intparams[m].mode);
540
541 fgets( ptr, MAXLINE, fp);
542 Tokenize( ptr, &words );
543
544 intparams[m].a = atof(words[1]) ;
545 intparams[m].p = atof(words[2]) ;
546 intparams[m].ksi = atof(words[3]) ;
547 intparams[m].q = atof(words[4]) ;
548 if (verbose) printf (" %s %f %f %f %f\n",words[0],
549 intparams[m].a,intparams[m].p,intparams[m].ksi,intparams[m].q);
550
551 // Ligne 6 - rayon de coupure potential SM
552
553 fgets( ptr, MAXLINE, fp);
554 Tokenize( ptr, &words );
555
556 intparams[m].dc1 = atof(words[1]) ;
557 intparams[m].dc2 = atof(words[2]) ;
558 intparams[m].r0 = atof(words[3]) ;
559
560
561 if (strcmp(intparams[m].mode,"metal") == 0) {
562 if (verbose) printf (" %s %f %f %f\n",words[0],
563 intparams[m].dc1,intparams[m].dc2,intparams[m].r0);
564 } else {
565 if (verbose) printf (" %s %f %f %f\n",words[0],
566 intparams[m].dc1,intparams[m].dc2,intparams[m].r0);
567 }
568
569
570 } else if (strcmp(intparams[m].typepot,"buck") == 0) {
571
572 if (verbose) printf(" %s %s %s %s\n",words[0],words[1],words[2],
573 intparams[m].typepot);
574
575 fgets( ptr, MAXLINE, fp);
576 Tokenize( ptr, &words );
577
578 intparams[m].abuck = atof(words[1]) ; intparams[m].rhobuck = atof(words[2]) ;
579 if (verbose) printf (" %s %f %f\n",words[0],intparams[m].abuck,intparams[m].rhobuck);
580
581 }
582
583 else if (strcmp(intparams[m].typepot,"buckPlusAttr") == 0) {
584
585 if (verbose) printf(" %s %s %s %s\n",words[0],words[1],words[2],
586 intparams[m].typepot);
587
588 fgets( ptr, MAXLINE, fp);
589 Tokenize( ptr, &words );
590
591 intparams[m].abuck = atof(words[1]) ; intparams[m].rhobuck = atof(words[2]) ;
592 if (verbose) printf (" %s %f %f\n",words[0],intparams[m].abuck,intparams[m].rhobuck);
593
594
595 fgets( ptr, MAXLINE, fp);
596 Tokenize( ptr, &words );
597
598 intparams[m].aOO = atof(words[1]) ; intparams[m].bOO = atof(words[2]) ;
599 intparams[m].r1OO = atof(words[3]) ;intparams[m].r2OO = atof(words[4]) ;
600 if (verbose) printf (" %s %f %f %f %f \n",words[0],intparams[m].aOO,
601 intparams[m].bOO,intparams[m].r1OO,intparams[m].r2OO);
602
603 }
604 if (verbose) printf (" intsm %d \n",intparams[m].intsm);
605
606 } // for maxintparam
607
608
609 /* ====================================================================
610 tables Parameters
611 ==================================================================== */
612
613 // Ligne 9 - rayon de coupure Electrostatique
614 if (test == 0) {
615 fgets(ptr,MAXLINE,fp);
616 if (verbose) printf ("%s\n",ptr);
617
618 fgets( ptr, MAXLINE, fp);
619 }
620 Tokenize( ptr, &words );
621
622 for (i=0 ; i<num_atom_types; i++) { params[i].cutsq = atof(words[1]); }
623 cutmax = atof(words[1]);
624 if (verbose) printf (" %s %f\n",words[0],params[0].cutsq);
625
626 // Ligne 9 - parametre pour les tableaux
627
628 fgets( ptr, MAXLINE, fp);
629 Tokenize( ptr, &words );
630
631 rmin = atof(words[1]) ; dr = atof(words[2]);
632 if (verbose) printf (" %s %f %f\n",words[0],rmin,dr);
633
634 kmax = int(cutmax*cutmax/(2.0*dr*rmin));
635 ds = cutmax*cutmax/static_cast<double>(kmax) ;
636 if (verbose) printf (" kmax %d et ds %f\n",kmax,ds);
637
638 /* ======================================================== */
639 fgets( ptr, MAXLINE, fp);
640 if (verbose) printf ("%s",ptr);
641
642 fgets( ptr, MAXLINE, fp);
643 Tokenize( ptr, &words );
644 Qstep = atoi(words[1]);
645 if (verbose) printf (" %s " BIGINT_FORMAT "\n",words[0],Qstep);
646
647 fgets( ptr, MAXLINE, fp);
648 Tokenize( ptr, &words );
649 loopmax = atoi(words[1]);
650 precision = atof(words[2]);
651 if (verbose) printf (" %s %d %f\n",words[0],loopmax,precision);
652
653 /* Param de coordination ============================================= */
654
655 fgets( ptr, MAXLINE, fp);
656 if (verbose) printf ("%s",ptr);
657
658 fgets( ptr, MAXLINE, fp);
659 Tokenize( ptr, &words );
660 r1Coord = atof(words[1]);
661 r2Coord = atof(words[2]);
662 if (verbose) printf (" %s %f %f\n",words[0],r1Coord,r2Coord);
663
664 /* Mode for QInit============================================= */
665 fgets( ptr, MAXLINE, fp);
666 if (verbose) printf ("%s",ptr);
667
668 fgets( ptr, MAXLINE, fp);
669 Tokenize( ptr, &words );
670 strcpy( QInitMode , words[1] );
671 if (strcmp(QInitMode,"true") == 0) QOxInit= atof(words[2]);
672 else QOxInit = 0.0;
673 if (verbose) printf (" %s %s %f\n",words[0],QInitMode,QOxInit);
674
675
676 /* Mode for QEq============================================= */
677
678 fgets( ptr, MAXLINE, fp);
679 if (verbose) printf ("%s",ptr);
680
681 fgets( ptr, MAXLINE, fp);
682 Tokenize( ptr, &words );
683 strcpy( QEqMode , words[1] );
684 if (verbose) printf (" %s %s\n",words[0],QEqMode);
685
686 fgets( ptr, MAXLINE, fp);
687
688 if (strcmp(QEqMode,"BulkFromSlab") == 0) {
689 Tokenize( ptr, &words );
690 zlim1QEq = atof(words[1]);
691 zlim2QEq = atof(words[2]);
692 if (verbose) printf (" %s %f %f\n",words[0],zlim1QEq,zlim2QEq);
693
694 } else if (strcmp(QEqMode,"Surface") == 0) {
695 Tokenize( ptr, &words );
696 zlim1QEq = atof(words[1]);
697 if (verbose) printf (" %s %f \n",words[0],zlim1QEq);
698
699 } else if (strcmp(QEqMode,"QEqAll") != 0 &&
700 strcmp(QEqMode,"QEqAllParallel") != 0 &&
701 strcmp(QEqMode,"Surface") != 0) {
702 error->all(FLERR,"The QEq Mode is not known. QEq mode should be :\n"
703 " Possible QEq modes | parameters\n"
704 " QEqAll | no parameters\n"
705 " QEqAllParallel | no parameters\n"
706 " Surface | zlim (QEq only for z>zlim)\n"
707 " BulkFromSlab | zlim1 zlim2 (QEq only for zlim1<z<zlim2)\n");
708 }
709
710 /* Bavard============================================= */
711
712 fgets( ptr, MAXLINE, fp);
713 if (verbose) printf ("%s",ptr);
714
715 fgets( ptr, MAXLINE, fp);
716 Tokenize( ptr, &words );
717 strcpy( Bavard , words[1] );
718 if (verbose) printf (" %s %s\n",words[0],Bavard);
719
720 // ---------------------------------------
721 // Writing the energy component.
722
723 fgets( ptr, MAXLINE, fp);
724 Tokenize( ptr, &words );
725 strcpy( writeenerg, words[1] );
726 if (strcmp (writeenerg,"true") == 0) { Nevery = atof(words[2]); }
727 else { Nevery = 0.0; }
728 if (verbose) printf (" %s %s %f\n",words[0],writeenerg,Nevery);
729
730 // ---------------------------------------
731 // Writing the chimical electronic potential.
732
733 fgets( ptr, MAXLINE, fp);
734 Tokenize( ptr, &words );
735 strcpy( writepot, words[1] );
736 if (strcmp (writepot,"true") == 0) { Neverypot = atof(words[2]); }
737 else { Neverypot = 0.0; }
738 if (verbose) printf (" %s %s %f\n",words[0],writepot,Neverypot);
739
740
741 /* ======================================================== */
742
743 /* deallocate helper storage */
744 for (i = 0; i < MAXTOKENS ; i++ ) free( words[i]);
745 free( words );
746 free( ptr );
747 fclose(fp);
748
749 // === Rayon de coupure premier voisins : 1,2*r0
750 for (i=0 ; i<num_atom_types ; i++) {
751 for (j=0 ; j<=i ; j++) {
752 m = intype[i][j];
753 if (m == 0) continue;
754 if (intparams[m].intsm == 0) continue;
755
756 intparams[m].neig_cut = 1.2*intparams[m].r0;
757 if (strcmp(intparams[m].typepot,"second_moment") == 0 )
758 if (verbose) printf (" Rc 1er voisin, typepot %s -> %f Ang\n",
759 intparams[m].typepot,intparams[m].neig_cut);
760 }
761 }
762
763 //A adapter au STO
764 for (i=1,ncov=params[0].sto*params[0].n0; i < nparams; ++i)
765 ncov = min(ncov,(params[1].sto)*(params[1].n0));
766
767 if (verbose) printf (" Parametre ncov = %f\n",ncov);
768 if (verbose) printf (" ********************************************* \n");
769 }
770
771 /* ----------------------------------------------------------------------
772 * COMPUTE
773 ---------------------------------------------------------------------- */
774
compute(int eflag,int vflag)775 void PairSMTBQ::compute(int eflag, int vflag)
776 {
777 int i,j,ii,jj,inum,jnum,m,gp;
778 tagint itag,jtag;
779 int itype,jtype;
780 int *ilist,*jlist,*numneigh,**firstneigh;
781
782 double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
783 double rsq,iq,jq,Eself;
784 double ecovtot,ErepOO,ErepMO,Eion,Ecoh;
785 double **tmp,**tmpAll,*nmol;
786 double dq,dqcov;
787
788 if (atom->nmax > nmax) {
789 memory->destroy(ecov);
790 memory->destroy(potmad);
791 memory->destroy(potself);
792 memory->destroy(potcov);
793 memory->destroy(sbcov);
794 memory->destroy(coord);
795 memory->destroy(sbmet);
796 memory->destroy(flag_QEq);
797 memory->destroy(qf);
798 memory->destroy(q1);
799 memory->destroy(q2);
800 memory->destroy(tab_comm);
801
802 nmax = atom->nmax;
803
804 memory->create(ecov,nmax,"pair:ecov");
805 memory->create(potmad,nmax,"pair:potmad");
806 memory->create(potself,nmax,"pair:potself");
807 memory->create(potcov,nmax,"pair:potcov");
808 memory->create(sbcov,nmax,"pair:sbcov");
809 memory->create(coord,nmax,"pair:coord");
810 memory->create(sbmet,nmax,"pair:sbmet");
811 memory->create(flag_QEq,nmax,"pair:flag_QEq");
812 memory->create(qf,nmax,"pair:qf");
813 memory->create(q1,nmax,"pair:q1");
814 memory->create(q2,nmax,"pair:q2");
815 memory->create(tab_comm,nmax,"pair:tab_comm");
816 }
817
818
819 evdwl = ecoul = ecovtot = ErepOO = ErepMO = Eion = 0.0;
820 Eself = 0.0;
821
822 ev_init(eflag,vflag);
823
824 double **x = atom->x;
825 double **f = atom->f;
826 double *q = atom->q;
827 tagint *tag = atom->tag;
828 int *type = atom->type;
829 int newton_pair = force->newton_pair;
830 int nlocal = atom->nlocal;
831
832 inum = list->inum;
833 ilist = list->ilist;
834 numneigh = list->numneigh;
835 firstneigh = list->firstneigh;
836
837 const bigint step = update->ntimestep;
838 if (step == 0 || ((Qstep > 0) && (step % Qstep == 0))) Charge();
839 // ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
840 // this is necessary to get sbcov or sbmet table in order to caclulate the covalent or metal bonding
841 if (Qstep == 0 || (step % Qstep != 0)) QForce_charge(0);
842
843
844 // Charges Communication
845 // ----------------------
846 forward(q) ; // reverse(q);
847
848 memory->create(nmol,nteam+1,"pair:nmol");
849 memory->create(tmp,nteam+1,7,"pair:tmp");
850 memory->create(tmpAll,nteam+1,7,"pair:tmpAll");
851
852
853 for (i=0; i<nteam+1; i++) {
854 nmol[i] = static_cast<double>(nQEqall[i]);
855 for (j=0; j<7; j++) { tmp[i][j] = 0.0; tmpAll[i][j] = 0.0; }
856 }
857
858
859 /* ------------------------------------------------------------------------
860 Energy component store in tmp[gp][:] with gp is # QEq group
861 0 -> ionic energy
862 1 -> coulombian energy
863 2 -> Electrosatic energy (ionic + Coulombian)
864 3 -> Short int. Ox-Ox
865 4 -> Short int. SMTB (repulsion)
866 5 -> Covalent energy SMTB
867 6 -> Somme des Q(i)²
868 ------------------------------------------------------------------------- */
869
870 /* -------------- N-body forces Calcul --------------- */
871
872 for (ii = 0; ii < inum; ii++) {
873 // ===============================
874 i = ilist[ii];
875 itag = tag[i];
876 itype = map[type[i]];
877 iq = q[i];
878 gp = flag_QEq[i];
879
880 xtmp = x[i][0];
881 ytmp = x[i][1];
882 ztmp = x[i][2];
883
884 // --- For atom i
885 tmp[gp][6] += iq*iq;
886
887
888 // self energy, only on i atom
889 // ---------------------------
890 Eself = self(¶ms[itype],iq);
891 tmp[gp][0] += Eself;
892 tmp[gp][2] += Eself;
893
894 if (evflag) ev_tally_full (i,0.0,2.0*Eself,0.0,0.0,0.0,0.0);
895
896
897 // N-body energy of i
898 // ---------------------
899 dq = fabs(params[itype].qform) - fabs(iq);
900 dqcov = dq*(2.0*ncov/params[itype].sto - dq);
901
902 ecov[i] = - sqrt(sbcov[i]*dqcov + sbmet[i]);
903 ecovtot += ecov[i];
904 tmp[gp][5] += ecov[i];
905
906 if (evflag) ev_tally_full(i,0.0,2.0*ecov[i],0.0,0.0,0.0,0.0);
907
908
909 // Coulombian Interaction
910 // -----------------------
911 evdwl = 2.0*Vself*iq*iq ;
912 tmp[gp][1] += Vself*iq*iq;
913 tmp[gp][2] += Vself*iq*iq;
914
915 if (evflag) ev_tally_full (i,0.0,evdwl,0.0,0.0,0.0,0.0);
916 evdwl = 0.0 ;
917
918
919 jlist = firstneigh[i];
920 jnum = numneigh[i];
921
922 for (jj = 0; jj < jnum; jj++) {
923 // ===============================
924 j = jlist[jj];
925 j &= NEIGHMASK;
926 jtype = map[type[j]];
927 jtag = tag[j]; jq = q[j];
928
929
930 // .......................................................................
931 if (itag > jtag) {
932 if ((itag+jtag) % 2 == 0) continue;
933 } else if (itag < jtag) {
934 if ((itag+jtag) % 2 == 1) continue;
935 } else {
936 if (x[j][2] < x[i][2]) continue;
937 if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
938 if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
939 }
940 // .......................................................................
941
942
943 // # of interaction
944 // ----------------
945 m = intype[itype][jtype];
946
947
948 delx = xtmp - x[j][0];
949 dely = ytmp - x[j][1];
950 delz = ztmp - x[j][2];
951 rsq = delx*delx + dely*dely + delz*delz;
952
953
954 // ---------------------------------
955 if (sqrt(rsq) > cutmax) continue;
956 // ---------------------------------
957
958
959 // Coulombian Energy
960 // ------------------
961 evdwl = 0.0 ; fpair = 0.0;
962 potqeq(i,j,iq,jq,rsq,fpair,eflag,evdwl);
963
964 tmp[gp][1] += evdwl;
965 tmp[gp][2] += evdwl;
966
967
968 // Coulombian Force
969 // -----------------
970 f[i][0] += delx*fpair;
971 f[i][1] += dely*fpair;
972 f[i][2] += delz*fpair;
973 f[j][0] -= delx*fpair;
974 f[j][1] -= dely*fpair;
975 f[j][2] -= delz*fpair;
976
977
978 if (evflag)
979 ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
980 evdwl = 0.0; fpair = 0.0 ;
981
982
983
984 // ---------------------
985 if (m == 0) continue;
986 // ---------------------
987
988 // ----------------------------------------------
989 if ( strcmp(intparams[m].typepot,"buck") == 0 ||
990 strcmp(intparams[m].typepot,"buckPlusAttr") ==0) {
991 // ----------------------------------------------
992
993 evdwl = 0.0; fpair =0.0;
994 rep_OO (&intparams[m],rsq,fpair,eflag,evdwl);
995 ErepOO += evdwl ;
996 tmp[gp][3] += evdwl;
997
998
999 // repulsion is pure two-body, sums up pair repulsive forces
1000 f[i][0] += delx*fpair;
1001 f[i][1] += dely*fpair;
1002 f[i][2] += delz*fpair;
1003 f[j][0] -= delx*fpair;
1004 f[j][1] -= dely*fpair;
1005 f[j][2] -= delz*fpair;
1006
1007
1008 if (evflag)
1009 ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
1010 evdwl = 0.0; fpair = 0.0 ;
1011
1012 } // ----------------------------------- Rep O-O
1013
1014 if (strcmp(intparams[m].typepot,"buckPlusAttr") == 0) {
1015 // ----------------------------------------------
1016
1017 evdwl = 0.0; fpair =0.0;
1018 Attr_OO (&intparams[m],rsq,fpair,eflag,evdwl);
1019 ErepOO += evdwl ;
1020 tmp[gp][3] += evdwl;
1021
1022
1023 // repulsion is pure two-body, sums up pair repulsive forces
1024 f[i][0] += delx*fpair;
1025 f[i][1] += dely*fpair;
1026 f[i][2] += delz*fpair;
1027 f[j][0] -= delx*fpair;
1028 f[j][1] -= dely*fpair;
1029 f[j][2] -= delz*fpair;
1030
1031
1032 if (evflag)
1033 ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
1034 evdwl = 0.0; fpair = 0.0 ;
1035
1036 } // ----------------------------------- Attr O-O
1037
1038
1039 // -----------------------------------------------------------------
1040 if (strcmp(intparams[m].typepot,"second_moment") != 0 ) continue;
1041 // -----------------------------------------------------------------
1042
1043
1044 if (sqrt(rsq) > intparams[m].dc2) continue;
1045 // -------------------------------------------
1046
1047 // Repulsion : Energy + force
1048 // ----------------------------
1049 evdwl = 0.0; fpair = 0.0 ;
1050 repulsive(&intparams[m],rsq,i,j,fpair,eflag,evdwl);
1051 ErepMO += evdwl;
1052 tmp[gp][4] += 2.0*evdwl;
1053
1054 f[i][0] += delx*fpair;
1055 f[i][1] += dely*fpair;
1056 f[i][2] += delz*fpair;
1057 f[j][0] -= delx*fpair;
1058 f[j][1] -= dely*fpair;
1059 f[j][2] -= delz*fpair;
1060
1061
1062 if (evflag)
1063 ev_tally(i,j,nlocal,newton_pair,2.0*evdwl,0.0,fpair,delx,dely,delz);
1064
1065 evdwl = 0.0 ; fpair = 0.0;
1066 // ----- ----- ----- ----- ----- -----
1067
1068 // Attraction : force
1069 // ------------------
1070 fpair = 0.0;
1071 f_att(&intparams[m], i, j, rsq, fpair) ;
1072
1073 f[i][0] += delx*fpair;
1074 f[i][1] += dely*fpair;
1075 f[i][2] += delz*fpair;
1076 f[j][0] -= delx*fpair;
1077 f[j][1] -= dely*fpair;
1078 f[j][2] -= delz*fpair;
1079
1080 if (evflag)
1081 ev_tally(i,j,nlocal,newton_pair,0.0,0.0,fpair,delx,dely,delz);
1082
1083
1084 } // --------------------------------- End j
1085
1086 } // ---------------------------------- End i
1087
1088
1089 if (vflag_fdotr) virial_fdotr_compute();
1090
1091
1092 for (i = 0; i < nteam+1; i++) {
1093 MPI_Allreduce(tmp[i],tmpAll[i],7,MPI_DOUBLE,MPI_SUM,world);
1094 }
1095
1096 // ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1097 if (me == 0 && fmod(static_cast<double>(step), Nevery) == 0.0 && strcmp(writeenerg,"true") == 0) {
1098
1099 ofstream fichierE;
1100
1101 if (step == 0) { fichierE.open ("Energy_component.txt", ios::out | ios::trunc) ;}
1102 else { fichierE.open ("Energy_component.txt", ios::out | ios::app) ;}
1103
1104 if (fichierE) fichierE<< setprecision(9) <<step;
1105
1106 for (gp = 0; gp < nteam+1; gp++) {
1107 if (nmol[gp] == 0) continue;
1108 if (fichierE) fichierE<< setprecision(9) <<" "<<gp<<" "<<nmol[gp]
1109 <<" "<<tmpAll[gp][2]<<" "<<tmpAll[gp][3]<<" "<<tmpAll[gp][4]+tmpAll[gp][5];
1110 }
1111 if (fichierE) fichierE<<endl;
1112 if (fichierE) fichierE.close();
1113 }
1114 // ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1115
1116 if (me == 0&& strcmp(Bavard,"false") != 0) {
1117 printf ("A la fin de Compute\n");
1118
1119 printf ("Nemton_pair : %d, evflag %d, tail_flag %d,vflag_fdotr %d\n",
1120 newton_pair,evflag,force->pair->tail_flag,vflag_fdotr);
1121 printf ("neighbor->includegroup %d\n",neighbor->includegroup);
1122
1123
1124
1125 for (gp=0; gp<nteam+1; gp++) { // -------------------------- Boucle sur les gp
1126
1127 printf ("Energie de la team %d -- %d atome -------(nteam = %d) \n ",gp,int(nmol[gp]),nteam);
1128
1129 if (nmol[gp] == 0) {
1130 printf (" ============================================ \n \n");
1131 continue;
1132 }
1133 printf ("Vself %f, Som q2 : %f, nmol %f\n",Vself,tmpAll[gp][6],nmol[gp]);
1134 // printf ("nmol %f\n",nmol[gp]);
1135 printf ("Energie coul tot : %f | %f par mol\n",tmpAll[gp][1],tmpAll[gp][1]/nmol[gp]);
1136 printf ("Energie ionique : %f | %f par mol\n",tmpAll[gp][0],tmpAll[gp][0]/nmol[gp]);
1137 printf ("Energie elect tot : %f | %f par mol\n",tmpAll[gp][2],tmpAll[gp][2]/nmol[gp]);
1138 printf ("Energie cp pair ox : %f | %f par mol\n",tmpAll[gp][3],tmpAll[gp][3]/nmol[gp]);
1139 printf ("Energie cp pair sm : %f | %f par mol\n",tmpAll[gp][4],tmpAll[gp][4]/nmol[gp]);
1140 printf ("Energie cov sm : %f | %f par mol\n",tmpAll[gp][5],tmpAll[gp][5]/nmol[gp]);
1141
1142 Ecoh = tmpAll[gp][2] + tmpAll[gp][3] + tmpAll[gp][4] + tmpAll[gp][5];
1143 printf ("Energie totale : %f | %f par mol\n",Ecoh,Ecoh/nmol[gp]);
1144 printf ("================================================= \n");
1145 printf (" \n");
1146
1147 } // ----------------------------------------------------- Boucle sur les gp
1148
1149
1150
1151
1152 } // ------------ Call me == 0
1153
1154 memory->destroy(nmol);
1155 memory->destroy(tmp);
1156 memory->destroy(tmpAll);
1157 }
1158
1159 /* ----------------------------------------------------------------------
1160 Partie Electrostatique
1161 ----------------------------------------------------------------------*/
1162
self(Param * param,double qi)1163 double PairSMTBQ::self(Param *param, double qi)
1164 {
1165 double self_tmp;
1166 double s1=param->chi, s2=param->dj;
1167
1168 self_tmp = qi*(s1+0.5*qi*s2);
1169
1170 return self_tmp;
1171 }
1172
1173 /* ---------------------------------------------------------------------- */
1174
qfo_self(Param * param,double qi)1175 double PairSMTBQ::qfo_self(Param *param, double qi)
1176 {
1177 double self_d;
1178 double s1 = param->chi;
1179 double s2 = param->dj;
1180
1181 self_d = 0.0 ;
1182 self_d = s1+qi*s2;
1183
1184 return self_d;
1185 }
1186
1187 /* ---------------------------------------------------------------------- */
1188 /* ---------------------------------------------------------------------- */
1189
tabqeq()1190 void PairSMTBQ::tabqeq()
1191 {
1192 int i,j,k,m,verbose;
1193 int nntype;
1194 double rc,s,r;
1195 double alf;
1196
1197 int ii;
1198 double za,zb,ra,rb,gam,dgam,dza,dzb,
1199 d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2,na,nb;
1200 double aCoeff,bCoeff,rcoupe,nang;
1201
1202 int n = atom->ntypes;
1203 int nlocal = atom->nlocal;
1204 int nghost = atom->nghost;
1205 nmax = atom->nmax;
1206
1207 verbose = 1;
1208 verbose = 0;
1209
1210 nntype = int((n+1)*n/2);
1211
1212 rc = cutmax ;
1213 alf = 0.3 ;
1214 // alf = 0.2 ;
1215
1216
1217 if (verbose) printf ("kmax %d, ds %f, nmax %d\n",kmax,ds,nmax);
1218 if (verbose) printf ("nlocal = %d, nghost = %d\n",nlocal,nghost);
1219 if (verbose) printf ("nntypes %d, kmax %d, rc %f, n %d\n",nntype,kmax,rc,n);
1220
1221 // allocate arrays
1222
1223 memory->create(coultype,n,n,"pair:intype");
1224 memory->create(potqn,kmax+5,"pair:potqn");
1225 memory->create(dpotqn,kmax+5,"pair:dpotqn");
1226 memory->create(fafb,kmax+5,nntype,"pair:fafb");
1227 memory->create(dfafb,kmax+5,nntype,"pair:dfafb");
1228 memory->create(fafbOxOxSurf,kmax+5,"pair:fafbOxOxSurf");
1229 memory->create(dfafbOxOxSurf,kmax+5,"pair:dfafbOxOxSurf");
1230 memory->create(fafbTiOxSurf,kmax+5,"pair:fafbTiOxSurf");
1231 memory->create(dfafbTiOxSurf,kmax+5,"pair:dfafbTiOxSurf");
1232
1233 memory->create(fafbOxOxBB,kmax+5,"pair:fafbOxOxBB");
1234 memory->create(dfafbOxOxBB,kmax+5,"pair:dfafbOxOxBB");
1235 memory->create(fafbTiOxBB,kmax+5,"pair:fafbTiOxB");
1236 memory->create(dfafbTiOxBB,kmax+5,"pair:dfafbTiOxBB");
1237
1238
1239 memory->create(ecov,nmax,"pair:ecov");
1240 memory->create(potmad,nmax,"pair:potmad");
1241 memory->create(potself,nmax,"pair:potself");
1242 memory->create(potcov,nmax,"pair:potcov");
1243 memory->create(sbcov,nmax,"pair:sbcov");
1244 memory->create(coord,nmax,"pair:coord");
1245 memory->create(sbmet,nmax,"pair:sbmet");
1246
1247 memory->create(flag_QEq,nmax,"pair:flag_QEq");
1248
1249 memory->create(qf,nmax,"pair:qf");
1250 memory->create(q1,nmax,"pair:q1");
1251 memory->create(q2,nmax,"pair:q2");
1252 memory->create(tab_comm,nmax,"pair:tab_comm");
1253
1254 memory->create(fct,31,"pair:fct");
1255
1256 // set interaction number: 0-0=0, 1-1=1, 0-1=1-0=2
1257
1258 m = 0; k = n;
1259 for (i = 0; i < n; i++) {
1260 for (j = 0; j < n; j++) {
1261 if (j == i) {
1262 coultype[i][j] = m;
1263 m += 1;
1264 } else if (j != i && j > i) {
1265 coultype[i][j] = k;
1266 k += 1;
1267 } else if (j != i && j < i) {
1268 coultype[i][j] = coultype[j][i];
1269 }
1270 if (verbose) printf ("i %d, j %d, coultype %d\n",i,j,coultype[i][j]);
1271 }
1272 }
1273
1274 // -------- Tabqn --------
1275
1276 // -------------------
1277 // Ouverture du fichier
1278 // ofstream fichier("tabqeq.txt", ios::out | ios::trunc) ;
1279 // -------------------
1280
1281 double mu;
1282
1283 mu = erfc(alf*rc)/rc ;
1284
1285 //if (fichier) fichier <<" r - potqn " <<endl ;
1286
1287 //-------------------------
1288 for (k=0; k < kmax+5; k++)
1289 //-------------------------
1290 {
1291 s = static_cast<double>(k)*ds ; r = sqrt(s);
1292 if (k==0) r=10e-30;
1293 potqn[k] = 14.4*(erfc(alf*r)/r - mu) ;
1294
1295 // $$$ Here is (1/r)*dE/dr
1296 dpotqn[k] = -14.4*( (erfc(alf*r)/(r*r) + 2.0*alf/MY_PIS/r*exp(-alf*alf*r*r))/r ) ;
1297 }
1298
1299 Vself = -14.4*(alf/MY_PIS + mu*0.5) ;
1300
1301 // --------------------
1302 // default arrays to zero
1303
1304 for (i = 0; i < kmax+5; i ++) {
1305 for (j = 0; j < nntype; j ++) {
1306 fafb[i][j] = 0.0;
1307 dfafb[i][j] = 0.0;
1308 }
1309 fafbOxOxSurf[i] = 0.0;
1310 fafbTiOxSurf[i] = 0.0;
1311 dfafbOxOxSurf[i] = 0.0;
1312 dfafbTiOxSurf[i] = 0.0;
1313
1314 fafbOxOxBB[i] = 0.0;
1315 fafbTiOxBB[i] = 0.0;
1316 dfafbOxOxBB[i] = 0.0;
1317 dfafbTiOxBB[i] = 0.0;
1318 }
1319
1320
1321 // Set Tabqeq
1322 double dij,ddij;
1323
1324 // -------------------
1325 // Ouverture du fichier
1326 //ofstream fichier("dtabqeq.txt", ios::out | ios::trunc) ;
1327 // -------------------
1328
1329 //if (fichier) fichier <<" k , r , fafb , dfafb , dfafb2 , dgam , d(1/r) , dpotqn" <<endl ;
1330
1331
1332 rcoupe = cutmax ;
1333 double cang ;
1334
1335 for (i = 0; i < n ; i++) {
1336 for (j = i; j < n ; j++) {
1337
1338 rc = cutmax; if (verbose) printf ("cutmax %f\n",cutmax);
1339 m = coultype[i][j] ;
1340 na = params[i].ne ;
1341 nb = params[j].ne ;
1342 za = params[i].dzeta ;
1343 zb = params[j].dzeta ;
1344 ra = params[i].R;
1345 rb = params[j].R;
1346
1347 ii = 0 ; nang =cang= 5.0 ;
1348 // --------------------------
1349 for (k = 0; k < kmax+5; k++)
1350 // --------------------------
1351 {
1352 gam = dgam = dza = dzb = d2zaa = d2zab =
1353 d2zbb = d2zra = d2zrb = d2gamr2 = 0.0 ;
1354 dij = 0.0 ;
1355
1356 s = static_cast<double>(k)*ds ; r = sqrt(s) ;
1357 if (k==0) r=10e-30;
1358
1359 gammas(na,nb,za,zb,r,gam,dgam,dza,dzb,
1360 d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
1361
1362 // --- Jij
1363
1364 dij = 14.4 * (1.0/r - static_cast<double>(gam));
1365 ddij = 14.4 * (-1.0/(r*r) - static_cast<double>(dgam)) ;
1366
1367 // Cutting Fonction
1368
1369 if (dij < 0.01 && ii==0)
1370 {
1371 ii=2;
1372 if (ii==2) if (verbose) printf ("rc : %f\n",r);
1373 rc = r ; ii=1 ;
1374 if ((rc+nang)>rcoupe) nang = rcoupe - rc ;
1375 bCoeff = (2*dij+ddij*nang)/(dij*nang);
1376 aCoeff = dij*exp(-bCoeff*rc) /square(nang);
1377 }
1378 if (r > rc) {dij = aCoeff *square(r- rc-nang) *exp(bCoeff*r);
1379 ddij = aCoeff*(r- rc-nang) *(2+bCoeff*(r-rc-nang))*exp(bCoeff*r);
1380 }
1381
1382 if (r > (rc+nang)) {dij = 0.0 ; ddij = 0.0;}
1383
1384 fafb[k][m] = potqn[k] - dij ;
1385 if (k == 1) fafb[0][m] = fafb[k][m] ;
1386
1387 dfafb[k][m] = dpotqn[k] - ddij/r ;
1388 }
1389
1390 // Make the table fafbOxOxSurf
1391 rc = cutmax;
1392 if (strcmp(params[i].nom,"O")==0 || strcmp(params[j].nom,"O")==0) {
1393 if (strcmp(params[i].nom,"O")==0) {
1394 ra = ROxSurf;
1395 za = (2.0*params[i].ne + 1.0)/(4.0*ra);}
1396 if (strcmp(params[j].nom,"O")==0) {
1397 rb = ROxSurf;
1398 zb = (2.0*params[j].ne + 1.0)/(4.0*rb); }
1399
1400 ii = 0 ; nang =cang= 5.0 ;
1401 // --------------------------
1402 for (k = 0; k < kmax+5; k++)
1403 // --------------------------
1404 {
1405 gam = dgam = dza = dzb = d2zaa = d2zab =
1406 d2zbb = d2zra = d2zrb = d2gamr2 = 0.0 ;
1407 dij = 0.0 ;
1408
1409 s = static_cast<double>(k)*ds ; r = sqrt(s) ;
1410 if (k==0) r=10e-30;
1411
1412 gammas(na,nb,za,zb,r,gam,dgam,dza,dzb,
1413 d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
1414
1415 // --- Jij
1416
1417 dij = 14.4 * (1.0/r - static_cast<double>(gam));
1418 ddij = 14.4 * (-1.0/(r*r) - static_cast<double>(dgam)) ;
1419
1420 if (dij < 0.01 && ii==0)
1421 {
1422 ii=2;
1423 if (ii==2) if (verbose) printf ("rc : %f\n",r);
1424 rc = r ; ii=1 ;
1425 if ((rc+nang)>rcoupe) nang = rcoupe - rc ;
1426 bCoeff = (2*dij+ddij*nang)/(dij*nang);
1427 aCoeff = dij*exp(-bCoeff*rc) /square(nang);
1428 }
1429 if (r > rc) {dij = aCoeff *square(r- rc-nang) *exp(bCoeff*r);
1430 ddij = aCoeff*(r- rc-nang) *(2+bCoeff*(r-rc-nang))*exp(bCoeff*r);
1431 }
1432
1433
1434 if (r > (rc+nang)) {dij = 0.0 ; ddij = 0.0;}
1435
1436 if (strcmp(params[i].nom,"O")==0 && strcmp(params[j].nom,"O")==0) {
1437 fafbOxOxSurf[k] = potqn[k] - dij ;
1438 if (k == 1) fafbOxOxSurf[0] = fafbOxOxSurf[k] ;
1439
1440 dfafbOxOxSurf[k] = dpotqn[k] - ddij/r ;
1441 }
1442 else { fafbTiOxSurf[k] = potqn[k] - dij ;
1443 if (k == 1) fafbTiOxSurf[0] = fafbTiOxSurf[k] ;
1444
1445 dfafbTiOxSurf[k] = dpotqn[k] - ddij/r ;}
1446
1447 }
1448
1449
1450 } //for k
1451 //end of make the table fafbOxOxSurf
1452
1453 // Makes the table fafbOxOxBB
1454 rc = cutmax;
1455 if (strcmp(params[i].nom,"O")==0 || strcmp(params[j].nom,"O")==0) {
1456 if (strcmp(params[i].nom,"O")==0) {
1457 ra = ROxBB;
1458 za = (2.0*params[i].ne + 1.0)/(4.0*ra);}
1459 if (strcmp(params[j].nom,"O")==0) {
1460 rb = ROxBB;
1461 zb = (2.0*params[j].ne + 1.0)/(4.0*rb); }
1462
1463
1464 ii = 0 ; nang =cang= 5.0 ;
1465 // --------------------------
1466 for (k = 0; k < kmax+5; k++)
1467 // --------------------------
1468 {
1469 gam = dgam = dza = dzb = d2zaa = d2zab =
1470 d2zbb = d2zra = d2zrb = d2gamr2 = 0.0 ;
1471 dij = 0.0 ;
1472
1473 s = static_cast<double>(k)*ds ; r = sqrt(s) ;
1474 if (k==0) r=10e-30;
1475
1476 gammas(na,nb,za,zb,r,gam,dgam,dza,dzb,
1477 d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
1478
1479 // --- Jij
1480
1481 dij = 14.4 * (1.0/r - static_cast<double>(gam));
1482 ddij = 14.4 * (-1.0/(r*r) - static_cast<double>(dgam)) ;
1483
1484 if (dij < 0.01 && ii==0) {
1485 ii=2;
1486 if (ii==2) if (verbose) printf ("rc : %f\n",r);
1487 rc = r ; ii=1 ;
1488 if ((rc+nang)>rcoupe) nang = rcoupe - rc ;
1489 bCoeff = (2*dij+ddij*nang)/(dij*nang);
1490 aCoeff = dij*exp(-bCoeff*rc) /square(nang);
1491 }
1492 if (r > rc) {
1493 dij = aCoeff *square(r- rc-nang) *exp(bCoeff*r);
1494 ddij = aCoeff*(r- rc-nang) *(2+bCoeff*(r-rc-nang))*exp(bCoeff*r);
1495 }
1496
1497
1498 if (r > (rc+nang)) {dij = 0.0 ; ddij = 0.0;}
1499
1500 if (strcmp(params[i].nom,"O")==0 && strcmp(params[j].nom,"O")==0) {
1501 fafbOxOxBB[k] = potqn[k] - dij ;
1502 if (k == 1) fafbOxOxBB[0] = fafbOxOxBB[k] ;
1503 dfafbOxOxBB[k] = dpotqn[k] - ddij/r ; }
1504 else { fafbTiOxBB[k] = potqn[k] - dij ;
1505 if (k == 1) fafbTiOxBB[0] = fafbTiOxBB[k] ;
1506 dfafbTiOxBB[k] = dpotqn[k] - ddij/r ;
1507 }
1508 }
1509
1510 } //for k
1511 //end of make the table fafbOxOxBB
1512
1513 }
1514 } //for i,j
1515
1516 //if (fichier) fichier.close() ;
1517
1518 }
1519
1520 /* ---------------------------------------------------------------------*/
1521
potqeq(int i,int j,double qi,double qj,double rsq,double & fforce,int,double & eng)1522 void PairSMTBQ::potqeq(int i, int j, double qi, double qj, double rsq,
1523 double &fforce, int /*eflag*/, double &eng)
1524 {
1525
1526 /* ===================================================================
1527 Coulombian energy calcul between i and j atoms
1528 with fafb table make in sm_table().
1529 fafb[i][j] : i is the table's step (r)
1530 j is the interaction's # (in intype[itype][jtype])
1531 dfafb is derivate energy (force)
1532 ==================================================================== */
1533
1534 int itype,jtype,l,m;
1535 double r,t1,t2,sds,xi,engBulk,engSurf,fforceBulk,fforceSurf,dcoordloc,dcoupureloc;
1536 double engBB,fforceBB, dIntfcoup2loc,iCoord,jCoord,iIntfCoup2,jIntfCoup2;
1537
1538 int *type = atom->type;
1539 // int n = atom->ntypes;
1540
1541 itype = map[type[i]];
1542 jtype = map[type[j]];
1543 m = coultype[itype][jtype];
1544
1545 r = rsq;
1546 sds = r/ds ; l = int(sds) ;
1547 xi = sds - static_cast<double>(l) ;
1548
1549
1550 iCoord=coord[i];
1551 jCoord=coord[j];
1552 iIntfCoup2= Intfcoup2(iCoord,coordOxBulk,0.15);
1553 jIntfCoup2= Intfcoup2(jCoord,coordOxBulk,0.15);
1554
1555 // ---- Energies Interpolation
1556
1557 t1 = fafb[l][m] + (fafb[l+1][m] - fafb[l][m])*xi;
1558 t2 = fafb[l+1][m] + (fafb[l+2][m] - fafb[l+1][m])*(xi-1.0);
1559
1560 engBulk = qi*qj*(t1 + (t2 - t1)*xi/2.0);
1561 eng=engBulk;
1562
1563
1564 // ---- Forces Interpolation
1565
1566 t1 = dfafb[l][m] + (dfafb[l+1][m] - dfafb[l][m])*xi;
1567 t2 = dfafb[l+1][m] + (dfafb[l+2][m] - dfafb[l+1][m])*(xi-1);
1568
1569
1570 fforce = - qi*qj*(t1 + (t2 - t1)*xi/2.0) ;
1571
1572
1573 if (strcmp(params[itype].nom,"O")==0 || strcmp(params[jtype].nom,"O")==0) {
1574
1575 if (strcmp(params[itype].nom,"O")==0 && strcmp(params[jtype].nom,"O")==0) {
1576 // between two oxygens
1577
1578 t1 = fafbOxOxSurf[l] + (fafbOxOxSurf[l+1] - fafbOxOxSurf[l])*xi;
1579 t2 = fafbOxOxSurf[l+1] + (fafbOxOxSurf[l+2] - fafbOxOxSurf[l+1])*(xi-1.0);
1580 engSurf = qi*qj*(t1 + (t2 - t1)*xi/2.0);
1581
1582 t1 = fafbOxOxBB[l] + (fafbOxOxBB[l+1] - fafbOxOxBB[l])*xi;
1583 t2 = fafbOxOxBB[l+1] + (fafbOxOxBB[l+2] - fafbOxOxBB[l+1])*(xi-1.0);
1584 engBB = qi*qj*(t1 + (t2 - t1)*xi/2.0);
1585
1586 eng= engBulk + (iCoord+jCoord-2*coordOxBulk)/(2*(coordOxBB-coordOxBulk)) *(engBB-engBulk)
1587 + (iIntfCoup2+jIntfCoup2)*((engBulk-engSurf)/(2*(coordOxBulk-coordOxSurf))
1588 - (engBB-engBulk)/(2*(coordOxBB-coordOxBulk))) ;
1589
1590
1591 // ---- Interpolation des forces
1592
1593 fforceBulk=fforce;
1594 t1 = dfafbOxOxSurf[l] + (dfafbOxOxSurf[l+1] - dfafbOxOxSurf[l])*xi;
1595 t2 = dfafbOxOxSurf[l+1] + (dfafbOxOxSurf[l+2] - dfafbOxOxSurf[l+1])*(xi-1);
1596 fforceSurf = - qi*qj*(t1 + (t2 - t1)*xi/2.0) ;
1597
1598 t1 = dfafbOxOxBB[l] + (dfafbOxOxBB[l+1] - dfafbOxOxBB[l])*xi;
1599 t2 = dfafbOxOxBB[l+1] + (dfafbOxOxBB[l+2] - dfafbOxOxBB[l+1])*(xi-1);
1600 fforceBB = - qi*qj*(t1 + (t2 - t1)*xi/2.0) ;
1601
1602 fforce= fforceBulk + (iCoord+jCoord-2*coordOxBulk)/(2*(coordOxBB-coordOxBulk))*(fforceBB-fforceBulk)
1603 + (iIntfCoup2+jIntfCoup2)*((fforceBulk-fforceSurf)/(2*(coordOxBulk-coordOxSurf))
1604 - (fforceBB-fforceBulk)/(2*(coordOxBB-coordOxBulk))) ;
1605
1606 } else { // between metal and oxygen
1607
1608 t1 = fafbTiOxSurf[l] + (fafbTiOxSurf[l+1] - fafbTiOxSurf[l])*xi;
1609 t2 = fafbTiOxSurf[l+1] + (fafbTiOxSurf[l+2] - fafbTiOxSurf[l+1])*(xi-1.0);
1610 engSurf = qi*qj*(t1 + (t2 - t1)*xi/2.0);
1611 t1 = fafbTiOxBB[l] + (fafbTiOxBB[l+1] - fafbTiOxBB[l])*xi;
1612 t2 = fafbTiOxBB[l+1] + (fafbTiOxBB[l+2] - fafbTiOxBB[l+1])*(xi-1.0);
1613 engBB = qi*qj*(t1 + (t2 - t1)*xi/2.0);
1614
1615 if (strcmp(params[jtype].nom,"O")==0) //the atom j is an oxygen
1616 { iIntfCoup2=jIntfCoup2;
1617 iCoord=jCoord; }
1618
1619 eng = engBulk + (engBulk-engSurf)/(coordOxBulk-coordOxSurf) * iIntfCoup2
1620 + (engBB-engBulk)/(coordOxBB-coordOxBulk) * (iCoord-coordOxBulk-iIntfCoup2);
1621
1622
1623 // ---- Forces Interpolation
1624
1625 fforceBulk=fforce;
1626 t1 = dfafbTiOxSurf[l] + (dfafbTiOxSurf[l+1] - dfafbTiOxSurf[l])*xi;
1627 t2 = dfafbTiOxSurf[l+1] + (dfafbTiOxSurf[l+2] - dfafbTiOxSurf[l+1])*(xi-1);
1628 fforceSurf = - qi*qj*(t1 + (t2 - t1)*xi/2.0) ;
1629
1630 t1 = dfafbTiOxBB[l] + (dfafbTiOxBB[l+1] - dfafbTiOxBB[l])*xi;
1631 t2 = dfafbTiOxBB[l+1] + (dfafbTiOxBB[l+2] - dfafbTiOxBB[l+1])*(xi-1);
1632 fforceBB = - qi*qj*(t1 + (t2 - t1)*xi/2.0) ;
1633
1634 dcoordloc = fcoupured(sqrt(r),r1Coord,r2Coord) ;
1635
1636
1637 dcoupureloc = fcoupured(iCoord,coordOxSurf,coordOxBulk) ;
1638 dIntfcoup2loc= fcoup2(iCoord,coordOxBulk,0.15)*dcoupureloc ;
1639 fforce = fforceBulk + 1/(coordOxBulk-coordOxSurf) * ((fforceBulk-fforceSurf)* iIntfCoup2
1640 - (engBulk-engSurf) *dIntfcoup2loc)
1641 + 1/(coordOxBB-coordOxBulk) * ((fforceBB-fforceBulk)*(iCoord-coordOxBulk- iIntfCoup2)
1642 - (engBB-engBulk) *(dcoordloc-dIntfcoup2loc));
1643
1644
1645 }
1646 }
1647 }
1648
1649 /* -------------------------------------------------------------------- */
1650
pot_ES(int i,int j,double rsq,double & eng)1651 void PairSMTBQ::pot_ES (int i, int j, double rsq, double &eng)
1652 {
1653
1654 /* ===================================================================
1655 Coulombian potential energy calcul between i and j atoms
1656 with fafb table make in sm_table().
1657 fafb[i][j] : i is the table's step (r)
1658 j is the interaction's # (in intype[itype][jtype])
1659 dfafb is derivate energy (force)
1660 ==================================================================== */
1661
1662 int itype,jtype,l,m;
1663 double r,t1,t2,sds,xi,engBulk,engSurf;
1664 double engBB,iCoord,jCoord,iIntfCoup2,jIntfCoup2;
1665
1666 int *type = atom->type;
1667 // int n = atom->ntypes;
1668
1669 itype = map[type[i]];
1670 jtype = map[type[j]];
1671 m = coultype[itype][jtype];
1672
1673 r = rsq;
1674 sds = r/ds ; l = int(sds) ;
1675 xi = sds - static_cast<double>(l) ;
1676
1677
1678 iCoord=coord[i];
1679 jCoord=coord[j];
1680 iIntfCoup2= Intfcoup2(iCoord,coordOxBulk,0.15);
1681 jIntfCoup2= Intfcoup2(jCoord,coordOxBulk,0.15);
1682
1683 // ---- Energies Interpolation
1684
1685 t1 = fafb[l][m] + (fafb[l+1][m] - fafb[l][m])*xi;
1686 t2 = fafb[l+1][m] + (fafb[l+2][m] - fafb[l+1][m])*(xi-1.0);
1687
1688
1689 eng = (t1 + (t2 - t1)*xi/2.0);
1690 engBulk=eng;
1691
1692
1693 if (itype==0 || jtype==0) {
1694
1695 if (itype==0 && jtype==0) { // between two oxygens
1696
1697 t1 = fafbOxOxSurf[l] + (fafbOxOxSurf[l+1] - fafbOxOxSurf[l])*xi;
1698 t2 = fafbOxOxSurf[l+1] + (fafbOxOxSurf[l+2] - fafbOxOxSurf[l+1])*(xi-1.0);
1699 engSurf = (t1 + (t2 - t1)*xi/2.0);
1700
1701 t1 = fafbOxOxBB[l] + (fafbOxOxBB[l+1] - fafbOxOxBB[l])*xi;
1702 t2 = fafbOxOxBB[l+1] + (fafbOxOxBB[l+2] - fafbOxOxBB[l+1])*(xi-1.0);
1703 engBB = (t1 + (t2 - t1)*xi/2.0);
1704
1705 eng= engBulk + (iCoord+jCoord-2*coordOxBulk)/(2*(coordOxBB-coordOxBulk))*(engBB-engBulk)
1706 + (iIntfCoup2+jIntfCoup2)*((engBulk-engSurf)/(2*(coordOxBulk-coordOxSurf))
1707 - (engBB-engBulk)/(2*(coordOxBB-coordOxBulk))) ;
1708
1709 } else { // between metal and oxygen
1710
1711 t1 = fafbTiOxSurf[l] + (fafbTiOxSurf[l+1] - fafbTiOxSurf[l])*xi;
1712 t2 = fafbTiOxSurf[l+1] + (fafbTiOxSurf[l+2] - fafbTiOxSurf[l+1])*(xi-1.0);
1713 engSurf = (t1 + (t2 - t1)*xi/2.0);
1714
1715 t1 = fafbTiOxBB[l] + (fafbTiOxBB[l+1] - fafbTiOxBB[l])*xi;
1716 t2 = fafbTiOxBB[l+1] + (fafbTiOxBB[l+2] - fafbTiOxBB[l+1])*(xi-1.0);
1717 engBB = (t1 + (t2 - t1)*xi/2.0);
1718
1719 if (jtype==0) { //the atom j is an oxygen
1720 iIntfCoup2=jIntfCoup2;
1721 iCoord=jCoord;
1722 }
1723
1724 eng = engBulk + (engBulk-engSurf)/(coordOxBulk-coordOxSurf)*iIntfCoup2
1725 + (engBB-engBulk)/(coordOxBB-coordOxBulk) * (iCoord-coordOxBulk-iIntfCoup2);
1726
1727
1728 }
1729
1730
1731
1732 }
1733
1734
1735 }
1736
1737 /* -------------------------------------------------------------------- */
1738
pot_ES2(int i,int j,double rsq,double & pot)1739 void PairSMTBQ::pot_ES2 (int i, int j, double rsq, double &pot)
1740 {
1741 int l,m,itype,jtype;
1742 double sds,xi,t1,t2,r;
1743
1744 int *type = atom->type;
1745
1746
1747 if (sqrt(rsq) > cutmax) return ;
1748
1749 itype = map[type[i]];
1750 jtype = map[type[j]];
1751 m = coultype[itype][jtype];
1752
1753 r = rsq ;
1754 sds = r/ds ; l = int(sds) ;
1755 xi = sds - static_cast<double>(l) ;
1756
1757 // ---- Energies Interpolation
1758
1759 t1 = fafb[l][m] + (fafb[l+1][m] - fafb[l][m])*xi;
1760 t2 = fafb[l+1][m] + (fafb[l+2][m] - fafb[l+1][m])*(xi-1.0);
1761
1762 pot = (t1 + (t2 - t1)*xi/2.0) ;
1763
1764 }
1765
1766 /* --------------------------------------------------------------------
1767 Oxygen-Oxygen Interaction
1768 -------------------------------------------------------------------- */
1769
rep_OO(Intparam * intparam,double rsq,double & fforce,int,double & eng)1770 void PairSMTBQ::rep_OO(Intparam *intparam, double rsq, double &fforce,
1771 int /*eflag*/, double &eng)
1772 {
1773 double r,tmp_exp,tmp;
1774 double A = intparam->abuck ;
1775 double rho = intparam->rhobuck ;
1776
1777 r = sqrt(rsq);
1778 tmp = - r/rho ;
1779 tmp_exp = exp( tmp );
1780
1781 eng = A * tmp_exp ;
1782
1783 fforce = A/rho * tmp_exp / r ; //( - )
1784
1785 }
1786
1787
Attr_OO(Intparam * intparam,double rsq,double & fforce,int,double & eng)1788 void PairSMTBQ::Attr_OO(Intparam *intparam, double rsq, double &fforce,
1789 int /*eflag*/, double &eng)
1790 {
1791 double r,tmp_exp;
1792 double aOO = intparam->aOO ;
1793 double bOO = intparam->bOO ;
1794 double r1OO = intparam->r1OO ;
1795 double r2OO = intparam->r2OO ;
1796
1797 r = sqrt(rsq);
1798 tmp_exp= exp( bOO* r);
1799 eng = aOO * tmp_exp* fcoupure(r,r1OO,r2OO);
1800
1801 fforce = - (aOO*bOO * tmp_exp * fcoupure(r,r1OO,r2OO)+ aOO*tmp_exp *fcoupured(r,r1OO,r2OO))/ r ; //( - )
1802
1803 }
1804
1805
1806 /* ----------------------------------------------------------------------
1807 covalente Interaction
1808 ----------------------------------------------------------------------*/
1809
1810
tabsm()1811 void PairSMTBQ::tabsm()
1812 {
1813 int k,m;
1814 double s,r,tmpb,tmpr,fcv,fcdv;
1815
1816 memory->create(tabsmb,kmax,maxintsm+1,"pair:tabsmb");
1817 memory->create(tabsmr,kmax,maxintsm+1,"pair:tabsmr");
1818 memory->create(dtabsmb,kmax,maxintsm+1,"pair:dtabsmb");
1819 memory->create(dtabsmr,kmax,maxintsm+1,"pair:dtabsmr");
1820
1821
1822 for (m = 0; m <= maxintparam; m++) {
1823
1824 if (intparams[m].intsm == 0) continue;
1825
1826 double rc1 = intparams[m].dc1;
1827 double rc2 = intparams[m].dc2;
1828 double A = intparams[m].a;
1829 double p = intparams[m].p;
1830 double Ksi = intparams[m].ksi;
1831 double q = intparams[m].q;
1832 double rzero = intparams[m].r0;
1833 int sm = intparams[m].intsm;
1834
1835
1836 for (k=0; k < kmax; k++)
1837 {
1838 s = static_cast<double>(k)*ds ; r = sqrt(s);
1839 if (k==0) r=10e-30;
1840 tmpb = exp( -2.0*q*(r/rzero - 1.0));
1841 tmpr = exp( -p*(r/rzero - 1.0));
1842
1843 if (r <= rc1)
1844 {
1845
1846 // -- Energy
1847 tabsmb[k][sm] = Ksi*Ksi * tmpb ;
1848 tabsmr[k][sm] = A * tmpr ;
1849
1850 // -- Force
1851 /* dtabsmb ne correspond pas vraiment a une force puisqu'il y a le /r
1852 (on a donc une unite force/distance). Le programme multiplie ensuite
1853 (dans le PairSMTBQ::compute ) dtabsmb par la projection du vecteur r
1854 sur un axe x (ou y ou z) pour determiner la composante de la force selon
1855 cette direction. Donc tout est ok au final. */
1856
1857 dtabsmb[k][sm] = - 2.0 *Ksi*Ksi* q/rzero * tmpb /r;
1858 dtabsmr[k][sm] = - A * p/rzero * tmpr/r ;
1859
1860 } // if
1861
1862 else if (r > rc1 && r <= rc2)
1863 {
1864
1865 // -- Energie
1866 fcv = fcoupure(r,intparams[sm].dc1,intparams[sm].dc2);
1867 tabsmb[k][sm] = fcv* Ksi*Ksi * tmpb ;
1868 tabsmr[k][sm] = fcv* A * tmpr ;
1869
1870 // -- Force
1871 /* dtabsmb ne correspond pas vraiment a une force puisqu'il y a le /r
1872 (on a donc une unite force/distance). Le programme multiplie ensuite
1873 (dans le PairSMTBQ::compute ) d tabsmb par la projection du vecteur
1874 r sur un axe x (ou y ou z) pour determiner la composante de la force
1875 selon cette direction. Donc tout est ok au final. */
1876
1877 fcdv = fcoupured(r,intparams[sm].dc1,intparams[sm].dc2);
1878 dtabsmb[k][sm] = (fcv*( - 2.0 *Ksi*Ksi* q/rzero * tmpb )+fcdv* Ksi*Ksi * tmpb )/r ;
1879 dtabsmr[k][sm] = (fcv*( - A * p/rzero * tmpr )+fcdv*A * tmpr )/r ;
1880
1881 }
1882
1883 else
1884 {
1885
1886 // -- Energie
1887 tabsmb[k][sm] = 0.0;
1888 tabsmr[k][sm] = 0.0;
1889
1890 // -- Force
1891 dtabsmb[k][sm] = 0.0;
1892 dtabsmr[k][sm] = 0.0;
1893
1894 }
1895
1896
1897
1898 } // for kmax
1899
1900
1901 } // for maxintparam
1902
1903 }
1904
1905
1906
1907
1908
1909 /* -------------------------------------------------------------- */
1910
repulsive(Intparam * intparam,double rsq,int,int,double & fforce,int,double & eng)1911 void PairSMTBQ::repulsive(Intparam *intparam, double rsq, int /*i*/, int /*j*/,
1912 double &fforce, int /*eflag*/, double &eng)
1913 {
1914
1915 /* ================================================
1916 rsq : square of ij distance
1917 fforce : repulsion force
1918 eng : repulsion energy
1919 eflag : Si oui ou non on calcule l'energie
1920 =================================================*/
1921
1922 int l;
1923 double r,sds,xi,t1,t2,dt1,dt2,sweet;
1924
1925 double rrcs = intparam->dc2;
1926 int sm = intparam->intsm;
1927
1928 // printf ("On rentre dans repulsive\n");
1929
1930
1931 r = rsq;
1932 if (sqrt(r) > rrcs) return ;
1933
1934 sds = r/ds ; l = int(sds) ;
1935 xi = sds - static_cast<double>(l) ;
1936
1937 t1 = tabsmr[l][sm] + (tabsmr[l+1][sm] - tabsmr[l][sm])*xi ;
1938 t2 = tabsmr[l+1][sm] + (tabsmr[l+2][sm] - tabsmr[l+1][sm])*(xi-1.0) ;
1939
1940 dt1 = dtabsmr[l][sm] + (dtabsmr[l+1][sm] - dtabsmr[l][sm])*xi ;
1941 dt2 = dtabsmr[l+1][sm] + (dtabsmr[l+2][sm] - dtabsmr[l+1][sm])*(xi-1.0) ;
1942
1943 if (strcmp(intparam->mode,"oxide") == 0)
1944 {
1945 fforce = - 2.0*(dt1 + (dt2 - dt1)*xi/2.0);
1946 eng = (t1 + (t2 - t1)*xi/2.0) ;
1947 }
1948 else if (strcmp(intparam->mode,"metal") == 0)
1949 {
1950 sweet = 1.0;
1951 fforce = - 2.0*(dt1 + (dt2 - dt1)*xi/2.0) * sweet ;
1952 eng = (t1 + (t2 - t1)*xi/2.0) * sweet ;
1953 }
1954
1955 }
1956
1957
1958 /* --------------------------------------------------------------------------------- */
1959
1960
attractive(Intparam * intparam,double rsq,int,int i,double,int,double)1961 void PairSMTBQ::attractive(Intparam *intparam, double rsq,
1962 int /*eflag*/, int i, double /*iq*/, int /*j*/, double /*jq*/)
1963 {
1964 int itype,l;
1965 double r,t1,t2,xi,sds;
1966 double sweet,mu;
1967
1968 double rrcs = intparam->dc2;
1969 int *type = atom->type;
1970 int sm = intparam->intsm;
1971
1972 itype = map[type[i]];
1973
1974 r = rsq;
1975 if (sqrt(r) > rrcs) return ;
1976
1977
1978 sds = r/ds ; l = int(sds) ;
1979 xi = sds - static_cast<double>(l) ;
1980
1981 t1 = tabsmb[l][sm] + (tabsmb[l+1][sm] - tabsmb[l][sm])*xi ;
1982 t2 = tabsmb[l+1][sm] + (tabsmb[l+2][sm] - tabsmb[l+1][sm])*(xi-1.0) ;
1983
1984
1985
1986 if (strcmp(intparam->mode,"oxide") == 0) {
1987 mu = 0.5*(sqrt(params[1].sto) + sqrt(params[0].sto));
1988
1989 // dq = fabs(params[itype].qform) - fabs(iq);
1990 // dqcov = dq*(2.0*ncov/params[itype].sto - dq);
1991
1992 sbcov[i] += (t1 + (t2 - t1)*xi/2.0) *params[itype].sto*mu*mu;
1993
1994 // if (i < 10) printf ("i %d, iq %f sbcov %f \n",i,iq,sbcov[i]);
1995
1996 if (sqrt(r)<r1Coord) { coord[i] += 1 ; }
1997 else if (sqrt(r)<r2Coord) { coord[i] += fcoupure(sqrt(r),r1Coord,r2Coord) ;}
1998
1999
2000 }
2001 else if (strcmp(intparam->mode,"metal") == 0) {
2002 sweet = 1.0;
2003 sbmet[i] += (t1 + (t2 - t1)*xi/2.0) * sweet ;
2004 }
2005
2006 }
2007
2008 /* ---------------------------------------------------------------------- */
2009
f_att(Intparam * intparam,int i,int j,double rsq,double & fforce)2010 void PairSMTBQ::f_att(Intparam *intparam, int i, int j,double rsq, double &fforce)
2011 {
2012 int itype,jtype,l;
2013 int *type = atom->type;
2014
2015 double r,sds,xi,dt1,dt2,dq,dqcovi,dqcovj;
2016 double fcov_ij,fcov_ji,sweet,iq,jq,mu;
2017
2018 int sm = intparam->intsm;
2019 double *q = atom->q;
2020
2021 itype = map[type[i]];
2022 jtype = map[type[j]];
2023 iq = q[i] ; jq = q[j];
2024
2025 r = rsq;
2026
2027 sds = r/ds ; l = int(sds) ;
2028 xi = sds - static_cast<double>(l) ;
2029
2030 dt1 = dtabsmb[l][sm] + (dtabsmb[l+1][sm] - dtabsmb[l][sm])*xi ;
2031 dt2 = dtabsmb[l+1][sm] + (dtabsmb[l+2][sm] - dtabsmb[l+1][sm])*(xi-1.0) ;
2032
2033 dq = fabs(params[itype].qform) - fabs(iq);
2034 dqcovi = dq*(2.0*ncov/params[itype].sto - dq);
2035
2036 dq = fabs(params[jtype].qform) - fabs(jq);
2037 dqcovj = dq*(2.0*ncov/params[jtype].sto - dq);
2038
2039 if (strcmp(intparam->mode,"oxide") == 0) {
2040 //------------------------------------------
2041 mu = 0.5*(sqrt(params[1].sto) + sqrt(params[0].sto));
2042 fcov_ij = (dt1 + (dt2 - dt1)*xi/2.0) * dqcovi *params[itype].sto*mu*mu;
2043 fcov_ji = (dt1 + (dt2 - dt1)*xi/2.0) * dqcovj *params[jtype].sto*mu*mu;
2044
2045 fforce = 0.5 * ( fcov_ij/sqrt(sbcov[i]*dqcovi + sbmet[i])
2046 + fcov_ji/sqrt(sbcov[j]*dqcovj + sbmet[j]) ) ;
2047 }
2048
2049 else if (strcmp(intparam->mode,"metal") == 0) {
2050 //-----------------------------------------------
2051 sweet = 1.0;
2052 fcov_ij = (dt1 + (dt2 - dt1)*xi/2.0) * sweet ;
2053
2054 fforce = 0.5 * fcov_ij*( 1.0/sqrt(sbcov[i]*dqcovi + sbmet[i])
2055 + 1.0/sqrt(sbcov[j]*dqcovj + sbmet[j]) ) ;
2056 }
2057
2058 }
2059
2060 /* ---------------------------------------------------------------------- */
2061
pot_COV(Param * param,int i,double & qforce)2062 void PairSMTBQ::pot_COV(Param *param, int i, double &qforce)
2063 {
2064 double iq,dq,DQ,sign;
2065
2066 double *q = atom->q;
2067 double qform = param->qform;
2068 double sto = param->sto;
2069
2070 sign = qform / fabs(qform);
2071 iq = q[i];
2072
2073 dq = fabs(qform) - fabs(iq);
2074 DQ = dq*(2.0*ncov/sto - dq);
2075
2076 if (fabs(iq) < 1.0e-7 || fabs(sbcov[i]) < 1.0e-7) {
2077 qforce = 0.0; }
2078 else {
2079 qforce = sign*sbcov[i]/sqrt(sbcov[i]*DQ + sbmet[i])*(ncov/sto - dq) ;
2080 }
2081
2082 }
2083
2084 /* ---------------------------------------------------------------------- */
2085
potmet(Intparam * intparam,double rsq,int i,double iq,int j,double jq)2086 double PairSMTBQ::potmet(Intparam *intparam, double rsq,
2087 int i, double iq, int j, double jq)
2088 {
2089 int l,itype,jtype;
2090 int *type = atom->type;
2091 double chi,sds,xi,t1,t2,r,dsweet,dq,dqcovi,dqcovj;
2092
2093 int sm = intparam->intsm;
2094 itype = map[type[i]];
2095 jtype = map[type[j]];
2096
2097 r = rsq;
2098 sds = r/ds ; l = int(sds) ;
2099 xi = sds - static_cast<double>(l) ;
2100
2101 t1 = tabsmb[l][sm] + (tabsmb[l+1][sm] - tabsmb[l][sm])*xi ;
2102 t2 = tabsmb[l+1][sm] + (tabsmb[l+2][sm] - tabsmb[l+1][sm])*(xi-1.0) ;
2103
2104 dq = fabs(params[itype].qform) - fabs(iq);
2105 dqcovi = dq*(2.0*ncov/params[itype].sto - dq);
2106
2107 dq = fabs(params[jtype].qform) - fabs(jq);
2108 dqcovj = dq*(2.0*ncov/params[jtype].sto - dq);
2109
2110 dsweet = 0.0;
2111 chi = (t1 + (t2 - t1)*xi/2.0) * dsweet *( 1.0/(2.0*sqrt(sbcov[i]*dqcovi+sbmet[i]))
2112 + 1.0/(2.0*sqrt(sbcov[j]*dqcovj+sbmet[j])) );
2113
2114 return chi;
2115 }
2116
2117
2118 /* ----------------------------------------------------------------------
2119 Cutting Function
2120 ----------------------------------------------------------------------*/
2121
2122
2123 /* -------------------------------------------------------------------- */
2124
fcoupure(double r,double rep_dc1,double rep_dc2)2125 double PairSMTBQ::fcoupure(double r, double rep_dc1, double rep_dc2)
2126 {
2127 double ddc,a3,a4,a5,x;
2128
2129 if (r<rep_dc1)
2130 {return 1;}
2131 else if (r> rep_dc2)
2132 {return 0;}
2133 else
2134 {
2135
2136 ddc = rep_dc2 - rep_dc1;
2137 x = r - rep_dc2;
2138
2139 a3 = -10/(ddc*ddc*ddc);
2140 a4 = -15/(ddc*ddc*ddc*ddc);
2141 a5 = -6/(ddc*ddc*ddc*ddc*ddc);
2142
2143 return x*x*x*(a3 + x*(a4 + x*a5));}
2144 }
2145
2146 /* ----------------------------------------------------------------------
2147 Derivate of cutting function
2148 ----------------------------------------------------------------------*/
2149
2150
2151 /* ----------------------------------------------------------------------- */
2152
fcoupured(double r,double rep_dc1,double rep_dc2)2153 double PairSMTBQ::fcoupured(double r, double rep_dc1, double rep_dc2)
2154 {
2155
2156 double ddc,a3,a4,a5,x;
2157
2158 if ( r>rep_dc1 && r<rep_dc2) {
2159 ddc = rep_dc2 - rep_dc1;
2160 x = r - rep_dc2;
2161
2162 a3 = -10/(ddc*ddc*ddc);
2163 a4 = -15/(ddc*ddc*ddc*ddc);
2164 a5 = -6/(ddc*ddc*ddc*ddc*ddc);
2165
2166 return x*x*(3*a3 + x*(4*a4 + 5*x*a5));}
2167 else {return 0;}
2168 }
2169
2170
2171 /* ----------------------------------------------------------------------
2172 cutting function for derive (force)
2173 ----------------------------------------------------------------------*/
2174
2175
2176 /* -------------------------------------------------------------------- */
2177
2178
2179
fcoup2(double c,double x,double delta)2180 double PairSMTBQ::fcoup2(double c, double x, double delta)
2181 {
2182 double dc;
2183
2184 if (c<x-delta)
2185 {return 1;}
2186 else if (c> x+delta)
2187 {return 0;}
2188 else
2189 {
2190 dc = c - x-delta;
2191 return dc*dc*(3*delta+dc)/(4*delta*delta*delta);
2192 }
2193 }
2194
2195 /* ----------------------------------------------------------------------
2196 Primitive of cutting function for derive (force)
2197 ----------------------------------------------------------------------*/
2198
2199
2200 /* -------------------------------------------------------------------- */
2201
Primfcoup2(double c,double x,double delta)2202 double PairSMTBQ::Primfcoup2(double c, double x, double delta)
2203 {
2204
2205 return (c*(c*c*c - 4* c*c* x - 4* (x - 2 *delta) * (x+delta)*(x+delta) +
2206 6* c *(x*x - delta*delta)))/(16* delta*delta*delta);
2207
2208 }
2209
2210
2211 /* ----------------------------------------------------------------------
2212 Integral of cutting function for derive (force) between x and c
2213 ----------------------------------------------------------------------*/
2214
2215
2216 /* -------------------------------------------------------------------- */
2217
Intfcoup2(double c,double x,double delta)2218 double PairSMTBQ::Intfcoup2(double c, double x, double delta)
2219 {
2220
2221 if (c<x-delta)
2222 {return c - x + delta + Primfcoup2(x-delta,x,delta) - Primfcoup2(x,x,delta) ;}
2223 else if (c> x+delta)
2224 {return Primfcoup2(x+delta,x,delta) - Primfcoup2(x,x,delta) ;}
2225 else
2226 {
2227 return Primfcoup2(c,x,delta) - Primfcoup2(x,x,delta) ;}
2228 }
2229
2230
2231 /* ---------------------------------------------------------------------
2232 Energy derivation respect charge Q
2233 --------------------------------------------------------------------- */
2234
QForce_charge(int loop)2235 void PairSMTBQ::QForce_charge(int loop)
2236 {
2237 int i,j,ii,jj,jnum;
2238 int itype,jtype,m,gp;
2239 double xtmp,ytmp,ztmp;
2240 double rsq;
2241 int *ilist,*jlist,*numneigh,**firstneigh;
2242 double iq,jq,fqi,fqj,fqij,fqij2,fqjj;
2243 int eflag = 0;
2244
2245
2246 double **x = atom->x;
2247 double *q = atom->q;
2248 int *type = atom->type;
2249 int step = update->ntimestep;
2250
2251 int inum = list->inum;
2252 ilist = list->ilist;
2253 numneigh = list->numneigh;
2254 firstneigh = list->firstneigh;
2255
2256 // loop over full neighbor list of my atoms
2257
2258 fqi = fqj = fqij = fqij2 = fqjj = 0.0;
2259
2260
2261 // ==================
2262 if (loop == 0) {
2263 // ==================
2264
2265 memset(sbcov,0,sizeof(double)*atom->nmax);
2266 memset(coord,0,sizeof(double)*atom->nmax);
2267 memset(sbmet,0,sizeof(double)*atom->nmax);
2268
2269 for (ii = 0; ii < inum; ii ++) {
2270 //--------------------------------
2271 i = ilist[ii];
2272 itype = map[type[i]];
2273
2274 gp = flag_QEq[i];
2275
2276 itype = map[type[i]];
2277 xtmp = x[i][0];
2278 ytmp = x[i][1];
2279 ztmp = x[i][2];
2280 iq = q[i];
2281
2282
2283
2284
2285 // two-body interactions
2286
2287 jlist = firstneigh[i];
2288 jnum = numneigh[i];
2289
2290 for (jj = 0; jj < jnum; jj++) {
2291 // -------------------------------
2292 j = jlist[jj];
2293 j &= NEIGHMASK;
2294
2295 jtype = map[type[j]];
2296 jq = q[j];
2297
2298 m = intype[itype][jtype];
2299 if (intparams[m].intsm == 0) continue ;
2300
2301
2302 const double delx = x[j][0] - xtmp;
2303 const double dely = x[j][1] - ytmp;
2304 const double delz = x[j][2] - ztmp;
2305 rsq = delx*delx + dely*dely + delz*delz;
2306
2307 // Covalente charge forces - sbcov initialization
2308 // ------------------------------------------------
2309 if (sqrt(rsq) > intparams[m].dc2) continue;
2310
2311 attractive (&intparams[m],rsq,eflag,i,iq,j,jq);
2312
2313
2314 } // ---------------------------------------------- for jj
2315
2316
2317 } // -------------------------------------------- ii
2318
2319 // ===============================================
2320 // Communicates the tables *sbcov and *sbmet
2321 // to calculate the N-body forces
2322 // ================================================
2323
2324 forward (sbcov) ; // reverse (sbcov);
2325 forward (coord) ; // reverse (coord);
2326 forward (sbmet) ; // reverse (sbmet);
2327
2328
2329 if (nteam == 0) return; //no oxide
2330 if (Qstep == 0 || (step % Qstep != 0)) return;
2331
2332 // =======================
2333 } // End of If(loop)
2334 // =======================
2335
2336
2337 // ===============================================
2338
2339 for (ii=0; ii<inum; ii++)
2340 {
2341 i = ilist[ii]; itype = map[type[i]];
2342 gp = flag_QEq[i]; iq = q[i];
2343
2344 potmad[i] = potself[i] = potcov[i] = 0.0 ;
2345
2346 if (gp == 0) continue;
2347
2348 xtmp = x[i][0];
2349 ytmp = x[i][1];
2350 ztmp = x[i][2];
2351
2352 fqi = 0.0 ;
2353
2354 // Madelung potential
2355 // --------------------
2356 potmad[i] += 2.0*Vself*iq ;
2357
2358 // charge force from self energy
2359 // -----------------------------
2360 fqi = qfo_self(¶ms[itype],iq);
2361 potself[i] = fqi ;
2362
2363
2364
2365 // Loop on Second moment neighbor
2366 // -------------------------------
2367
2368 jlist = firstneigh[i];
2369 jnum = numneigh[i];
2370
2371 for (jj = 0; jj < jnum; jj++)
2372 {
2373 j = jlist[jj];
2374 j &= NEIGHMASK;
2375 jtype = map[type[j]];
2376 m = intype[itype][jtype];
2377 jq = q[j];
2378
2379
2380 const double delx = x[j][0] - xtmp;
2381 const double dely = x[j][1] - ytmp;
2382 const double delz = x[j][2] - ztmp;
2383 rsq = delx*delx + dely*dely + delz*delz;
2384
2385 // long range q-dependent
2386 if (sqrt(rsq) > cutmax) continue;
2387
2388 // 1/r charge forces
2389 // --------------------
2390 fqij = 0.0;
2391 // pot_ES2 (i,j,rsq,fqij2);
2392 pot_ES (i,j,rsq,fqij);
2393
2394 potmad[i] += jq*fqij ;
2395
2396
2397 } // ------ jj
2398
2399 fqi = 0.0;
2400 pot_COV (¶ms[itype],i,fqi) ;
2401 potcov[i] = fqi ;
2402
2403
2404 } // ------- ii
2405
2406
2407 }
2408
2409 /* ---------------------------------------------------------------------- */
2410
Charge()2411 void PairSMTBQ::Charge()
2412 {
2413 int i,ii,iloop,itype,gp,m;
2414 int *ilist;
2415 double heatpq,qmass,dtq,dtq2,qtot,qtotll;
2416 double t_init,t_end,dt;
2417
2418 double *Transf,*TransfAll;
2419
2420 double *q = atom->q;
2421 double **x = atom->x;
2422 int *type = atom->type;
2423 int step = update->ntimestep;
2424
2425 int inum = list->inum;
2426 ilist = list->ilist;
2427
2428
2429 if (me == 0) t_init = MPI_Wtime();
2430 if (step == 0) cluster = 0;
2431
2432 // ---------------------------
2433 // Mise en place des groupes
2434 // ---------------------------
2435
2436 if (strcmp(QEqMode,"BulkFromSlab") == 0)
2437 groupBulkFromSlab_QEq();
2438 else if (strcmp(QEqMode,"QEqAll") == 0)
2439 groupQEqAll_QEq();
2440 else if (strcmp(QEqMode,"QEqAllParallel") == 0)
2441 groupQEqAllParallel_QEq();
2442 else if (strcmp(QEqMode,"Surface") == 0)
2443 groupSurface_QEq();
2444
2445
2446
2447 if (nteam+1 != cluster) {
2448 memory->destroy(nQEqall);
2449 memory->destroy(nQEqaall);
2450 memory->destroy(nQEqcall);
2451
2452 cluster = nteam+1;
2453 memory->create(nQEqall,nteam+1,"pair:nQEq");
2454 memory->create(nQEqaall,nteam+1,"pair:nQEqa");
2455 memory->create(nQEqcall,nteam+1,"pair:nQEqc");
2456 }
2457 // ---------------------------
2458
2459
2460 std::vector<double> enegtotall(nteam+1),enegchkall(nteam+1),enegmaxall(nteam+1);
2461 std::vector<double> qtotcll(nteam+1),qtotall(nteam+1),qtota(nteam+1),qtotc(nteam+1);
2462 std::vector<double> sigmaa(nteam+1),sigmac(nteam+1),sigmaall(nteam+1),sigmacll(nteam+1);
2463 std::vector<int> end(nteam+1), nQEq(nteam+1),nQEqc(nteam+1),nQEqa(nteam+1);
2464
2465
2466 iloop = 0;
2467
2468
2469 heatpq = 0.07;
2470 qmass = 0.000548580;
2471 dtq = 0.0006; // 0.0006
2472 dtq2 = 0.5*dtq*dtq/qmass;
2473
2474 std::vector<double> enegchk(nteam+1),enegtot(nteam+1),enegmax(nteam+1);
2475
2476
2477
2478 for (i=0; i<nteam+1; i++) {
2479 nQEq[i] = nQEqa[i] = nQEqc[i] = 0;
2480 nQEqall[i] = nQEqcall[i] = nQEqaall[i] = end[i]= 0;
2481 enegchk[i] = enegtot[i] = enegmax[i] = 0.0;
2482 qtota[i] = qtotc[i] = qtotall[i] = qtotcll[i] = 0.0;
2483 sigmaall[i] = sigmacll[i] = 0.0;
2484 }
2485 qtot = qtotll = 0.0 ;
2486
2487
2488 for (ii = 0; ii < inum; ii++) {
2489 i = ilist[ii]; itype = map[type[i]];
2490 gp = flag_QEq[i];
2491 q1[i] = q2[i] = qf[i] = 0.0;
2492
2493 qtot += q[i] ;
2494 nQEq[gp] += 1;
2495 if (itype != 0) { qtotc[gp] += q[i]; nQEqc[gp] += 1; }
2496 if (itype == 0) { qtota[gp] += q[i]; nQEqa[gp] += 1; }
2497 }
2498
2499 MPI_Allreduce(nQEq.data(),nQEqall,nteam+1,MPI_INT,MPI_SUM,world);
2500 MPI_Allreduce(nQEqc.data(),nQEqcall,nteam+1,MPI_INT,MPI_SUM,world);
2501 MPI_Allreduce(nQEqa.data(),nQEqaall,nteam+1,MPI_INT,MPI_SUM,world);
2502
2503 MPI_Allreduce(qtotc.data(),qtotcll.data(),nteam+1,MPI_DOUBLE,MPI_SUM,world);
2504 MPI_Allreduce(qtota.data(),qtotall.data(),nteam+1,MPI_DOUBLE,MPI_SUM,world);
2505 MPI_Allreduce(&qtot,&qtotll,1,MPI_DOUBLE,MPI_SUM,world);
2506
2507
2508 // ---------------------------------
2509 // Init_charge(nQEq,nQEqa,nQEqc);
2510 // ---------------------------------
2511
2512 if (update->ntimestep == 0 && (strcmp(QInitMode,"true") == 0)) {
2513 //Carefull here it won't be fine if there are more than 2 species!!!
2514 QOxInit=max(QOxInit, -0.98* params[1].qform *nQEqcall[gp]/nQEqaall[gp]) ;
2515
2516 for (ii = 0; ii < inum; ii++) {
2517 i = ilist[ii]; itype = map[type[i]];
2518 gp = flag_QEq[i];
2519
2520 // if (gp == 0) continue;
2521
2522 if (itype == 0) q[i] = QOxInit ;
2523 if (itype > 0) q[i] = -QOxInit * nQEqaall[gp]/nQEqcall[gp];
2524 }
2525 }
2526
2527 if (nteam == 0 || Qstep == 0) return;
2528 if (step % Qstep != 0) return;
2529 // --------------------------------------
2530
2531 // ----
2532
2533 // ----
2534
2535
2536 if (me == 0 && strcmp(Bavard,"false") != 0) {
2537 for (gp = 0; gp < nteam+1; gp++) {
2538 printf (" ++++ Groupe %d - Nox %d Ncat %d\n",gp,nQEqaall[gp],nQEqcall[gp]);
2539 if (nQEqcall[gp] !=0 && nQEqaall[gp] !=0 )
2540 printf (" neutralite des charges %f\n qtotc %f qtota %f\n",
2541 qtotll,qtotcll[gp]/nQEqcall[gp],qtotall[gp]/nQEqaall[gp]);
2542 printf (" ---------------------------- \n");}
2543 }
2544
2545 // ======================= Tab transfert ==================
2546 // Transf[gp] = enegtot[gp]
2547 // Transf[gp+cluster] = Qtotc[gp]
2548 // Transf[gp+2cluster] = Qtota[gp]
2549 // Transf[3cluster] = Qtot
2550 // -------------------------------------------------------
2551 memory->create(Transf,3*cluster+1,"pair:Tranf");
2552 memory->create(TransfAll,3*cluster+1,"pair:TranfAll");
2553 // ========================================================
2554
2555
2556 // --------------------------------------------
2557 for (iloop = 0; iloop < loopmax; iloop ++) {
2558 // --------------------------------------------
2559
2560 qtot = qtotll = Transf[3*cluster] = 0.0 ;
2561 for (gp=0; gp<nteam+1; gp++) {
2562 Transf[gp] = Transf[gp+cluster] = Transf[gp+2*cluster] = 0.0;
2563 TransfAll[gp] = TransfAll[gp+cluster] = TransfAll[gp+2*cluster] = 0.0;
2564 enegtot[gp] = enegtotall[gp] = enegchkall[gp] = enegmaxall[gp] = 0.0 ;
2565 }
2566
2567 for (ii = 0; ii < inum; ii++) {
2568 i = ilist[ii];
2569 q1[i] += qf[i]*dtq2 - heatpq*q1[i];
2570 q[i] += q1[i];
2571
2572 Transf[3*cluster] += q[i];
2573 itype = map[type[i]];
2574 gp = flag_QEq[i];
2575 if (itype == 0) Transf[gp+2*cluster] += q[i];
2576 if (itype != 0) Transf[gp+cluster] += q[i];
2577 }
2578
2579 // Communication des charges
2580 // ---------------------------
2581 forward(q) ; // reverse(q);
2582
2583
2584 // Calcul des potential
2585 // ----------------------
2586 QForce_charge(iloop);
2587
2588 // exit(1);
2589
2590 for (ii = 0; ii < inum; ii++)
2591 {
2592 i = ilist[ii]; itype = map[type[i]];
2593 gp = flag_QEq[i];
2594
2595 qf[i] = 0.0;
2596 // AK: chimet is not set anywhere
2597 qf[i] = potself[i]+potmad[i]+potcov[i]; // +chimet[i];
2598 Transf[gp] += qf[i];
2599 }
2600
2601 MPI_Allreduce(Transf,TransfAll,3*cluster+1,MPI_DOUBLE,MPI_SUM,world);
2602
2603 for (i = 0; i < nteam+1; i++) {
2604
2605 if (nQEqall[i] !=0) TransfAll[i] /= static_cast<double>(nQEqall[i]);
2606 enegchk[i] = enegmax[i] = 0.0;
2607 }
2608
2609 for (ii = 0; ii < inum; ii++) {
2610 i = ilist[ii];
2611 itype = map[type[i]];
2612 gp = flag_QEq[i];
2613
2614 if (gp == 0) continue;
2615
2616 q2[i] = TransfAll[gp] - qf[i];
2617
2618 enegmax[gp] = MAX(enegmax[gp],fabs(q2[i]));
2619 enegchk[gp] += fabs(q2[i]);
2620 qf[i] = q2[i];
2621
2622 } // Boucle local
2623
2624 MPI_Allreduce(enegchk.data(),enegchkall.data(),nteam+1,MPI_DOUBLE,MPI_SUM,world);
2625 MPI_Allreduce(enegmax.data(),enegmaxall.data(),nteam+1,MPI_DOUBLE,MPI_MAX,world);
2626
2627
2628 for (gp = 0; gp < nteam+1; gp++) {
2629 if (nQEqall[gp] !=0) {
2630 enegchk[gp] = enegchkall[gp]/static_cast<double>(nQEqall[gp]);
2631 enegmax[gp] = enegmaxall[gp];
2632 }
2633 }
2634
2635 // -----------------------------------------------------
2636 // Convergence Test
2637 // -----------------------------------------------------
2638
2639 m = 0;
2640 for (gp = 1; gp < nteam+1; gp++) {
2641 if (enegchk[gp] <= precision && enegmax[gp] <= 100.0*precision) end[gp] = 1;
2642 }
2643 for (gp = 1; gp < nteam+1; gp++) { m += end[gp] ; }
2644
2645 if (m == nteam) { break; }
2646 // -----------------------------------------------------
2647 // -----------------------------------------------------
2648
2649 for (ii = 0; ii < inum; ii++) {
2650 i = ilist[ii];
2651 q1[i] += qf[i]*dtq2 - heatpq*q1[i];
2652 }
2653
2654 // --------------------------------------------
2655 } // -------------------------------- iloop
2656 // --------------------------------------------
2657
2658
2659 // =======================================
2660 // Charge Communication.
2661 // =======================================
2662 forward(q); // reverse(q);
2663
2664 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
2665 // ==========================================
2666 // Ecriture des potentials dans un fichier
2667 // ==========================================
2668
2669 if (strcmp(writepot,"true") == 0 && fmod(static_cast<double>(step), Neverypot) == 0.0) {
2670
2671 ofstream fichierpot("Electroneg_component.txt", ios::out | ios::trunc) ;
2672
2673 for (ii = 0; ii < inum; ii++) {
2674 i = ilist[ii];
2675 itype = map[type[i]];
2676 gp = flag_QEq[i];
2677
2678
2679 if (fichierpot) fichierpot<< setprecision(9) <<i <<" "<<itype<<" "<<x[i][0]<<" "<<x[i][1]
2680 <<" "<<x[i][2]<<" "<<q[i]<<" "<<potself[i] + potmad[i]<<" "<<potcov[i]
2681 <<" "<<sbcov[i]<<" "<<TransfAll[gp]<<endl;
2682
2683 }
2684 if (fichierpot) fichierpot.close() ;
2685 }
2686
2687 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
2688
2689
2690 // Statistique (ecart type)
2691 // ------------------------
2692 for (i=0; i<nteam+1; i++) {
2693 if (nQEqcall[i] !=0)
2694 { TransfAll[i+cluster] /= static_cast<double>(nQEqcall[i]) ;
2695 TransfAll[i+2*cluster] /= static_cast<double>(nQEqaall[i]) ;}
2696 sigmaa[i] = sigmac[i] = 0.0;
2697 }
2698
2699 qtot = 0.0 ;
2700 for (ii = 0; ii < inum; ii++)
2701 {
2702 i = ilist[ii];
2703 itype = map[type[i]];
2704 gp = flag_QEq[i];
2705 // qtot += q[i];
2706
2707 if (gp == 0) continue;
2708
2709 if (itype == 0) sigmaa[gp] += (q[i]-TransfAll[gp+2*cluster])*(q[i]-TransfAll[gp+2*cluster]);
2710 if (itype == 1) sigmac[gp] += (q[i]-TransfAll[gp+cluster])*(q[i]-TransfAll[gp+cluster]);
2711 }
2712
2713 MPI_Allreduce(sigmaa.data(),sigmaall.data(),nteam+1,MPI_DOUBLE,MPI_SUM,world);
2714 MPI_Allreduce(sigmac.data(),sigmacll.data(),nteam+1,MPI_DOUBLE,MPI_SUM,world);
2715
2716 for (gp = 1; gp < nteam+1; gp++) {
2717 sigmaall[gp] = sqrt(sigmaall[gp]/static_cast<double>(nQEqaall[gp])) ;
2718 sigmacll[gp] = sqrt(sigmacll[gp]/static_cast<double>(nQEqcall[gp])) ;
2719 }
2720
2721
2722
2723 if (me == 0 && strcmp(Bavard,"false") != 0) {
2724 for (gp = 0; gp < nteam+1; gp++) {
2725 printf (" -------------- Groupe %d -----------------\n",gp);
2726 printf (" qtotc %f(+- %f) qtota %f(+- %f)\n",
2727 TransfAll[gp+cluster],sigmacll[gp],TransfAll[gp+2*cluster],sigmaall[gp]);
2728 printf (" Potentiel elec total : %f\n iloop %d, qtot %f\n",TransfAll[gp],iloop,TransfAll[3*cluster]);
2729 printf (" convergence : %f - %f\n",enegchk[gp],enegmax[gp]);
2730 }
2731
2732 t_end = MPI_Wtime();
2733 dt = t_end - t_init;
2734 printf (" temps dans charges : %f seconde. \n",dt);
2735 printf (" ======================================================== \n");
2736 }
2737
2738 // ============== Destroy Tab
2739 memory->destroy(Transf);
2740 memory->destroy(TransfAll);
2741
2742 }
2743
2744 /* ---------------------------------------------------------------------- */
2745
groupBulkFromSlab_QEq()2746 void PairSMTBQ::groupBulkFromSlab_QEq()
2747 { int ii,i;
2748 double **x=atom->x;
2749 int *ilist;
2750 double ztmp;
2751 int inum = list->inum;
2752 ilist = list->ilist;
2753
2754 for (ii = 0; ii < inum; ii++)
2755 {
2756 i = ilist[ii];
2757 ztmp = x[i][2];
2758 if (ztmp>zlim1QEq && ztmp< zlim2QEq)
2759 flag_QEq[i]=1;
2760 else
2761 flag_QEq[i]=0;
2762
2763 nteam=1;
2764
2765 }
2766
2767 }
2768
2769 // ----------------------------------------------
2770
groupQEqAll_QEq()2771 void PairSMTBQ::groupQEqAll_QEq()
2772 { int ii,i;
2773 int *ilist;
2774 int inum = list->inum;
2775 ilist = list->ilist;
2776
2777 nteam=1;
2778
2779 for (ii = 0; ii < inum; ii++)
2780 {
2781 i= ilist[ii];
2782 flag_QEq[i]=1;
2783 }
2784
2785 }
2786
2787 // ----------------------------------------------
2788
groupSurface_QEq()2789 void PairSMTBQ::groupSurface_QEq()
2790 { int ii,i;
2791 double **x=atom->x;
2792 int *ilist;
2793 double ztmp;
2794 int inum = list->inum;
2795 ilist = list->ilist;
2796
2797 for (ii = 0; ii < inum; ii++)
2798 {
2799 i = ilist[ii];
2800 ztmp = x[i][2];
2801 if (ztmp>zlim1QEq)
2802 flag_QEq[i]=1;
2803 else
2804 flag_QEq[i]=0;
2805
2806 nteam=1;
2807
2808 }
2809
2810 }
2811
2812
groupQEqAllParallel_QEq()2813 void PairSMTBQ::groupQEqAllParallel_QEq()
2814 {
2815 int ii,i,jj,j,kk,k,itype,jtype,ktype,jnum,m,gp,zz,z,kgp;
2816 int iproc; // ,team_elt[10][nproc],team_QEq[10][nproc][5];
2817 int **team_elt,***team_QEq;
2818 int *ilist,*jlist,*numneigh,**firstneigh,ngp,igp;
2819 double delr[3],xtmp,ytmp,ztmp,rsq;
2820 int **flag_gp, *nelt, **tab_gp;
2821 int QEq,*QEqall;
2822
2823 double **x = atom->x;
2824 int *type = atom->type;
2825 const int nlocal = atom->nlocal;
2826 const int nghost = atom->nghost;
2827 const int nall = nlocal + nghost;
2828 int inum = list->inum;
2829 ilist = list->ilist;
2830 numneigh = list->numneigh;
2831 firstneigh = list->firstneigh;
2832
2833
2834 // +++++++++++++++++++++++++++++++++++++++++++++++++
2835 // On declare et initialise nos p'tits tableaux
2836 // +++++++++++++++++++++++++++++++++++++++++++++++++
2837
2838 int **tabtemp,**Alltabtemp, *gptmp, *Allgptmp;
2839
2840 memory->create(tabtemp,10*nproc+10,nproc,"pair:tabtemp");
2841 memory->create(Alltabtemp,10*nproc+10,nproc,"pair:Alltabtemp");
2842 memory->create(gptmp,10*nproc+10,"pair:gptmp");
2843 memory->create(Allgptmp,10*nproc+10,"pair:Allgptmp");
2844
2845 memory->create(flag_gp,nproc,nall,"pair:flag_gp");
2846 memory->create(nelt,nall,"pair:nelt");
2847 memory->create(tab_gp,10,nall,"pair:flag_gp");
2848
2849 memory->create(team_elt,10,nproc,"pair:team_elt");
2850 memory->create(team_QEq,10,nproc,5,"pair:team_QEq");
2851 memory->create(QEqall,nproc,"pair:QEqall");
2852
2853 for (i = 0; i < nall ; i++) { flag_QEq[i] = 0; }
2854 for (i = 0; i < 10*nproc; i++) {
2855 gptmp[i] = 0; Allgptmp[i] = 0;
2856 for (j=0;j<nproc;j++) { tabtemp[i][j] = 0;
2857 Alltabtemp[i][j] = 0;}
2858 }
2859 for (i = 0; i < 10; i++) {
2860 for (k = 0; k < nall; k++) { tab_gp[i][k] = 0;
2861 if (i == 0) nelt[k] = 0;
2862 }
2863 for (j = 0; j < nproc; j++) {
2864 team_elt[i][j] = 0;
2865 for (k = 0; k < 5; k++) { team_QEq[i][j][k] = 0; }
2866 }
2867 }
2868
2869 QEq = 0;
2870
2871
2872 // printf ("groupeQEq me %d - nloc %d nghost %d boite %d\n",
2873 // me,nlocal,nghost,nall);
2874
2875 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++
2876 // On identifie les atomes rentrant dans le schema QEq +
2877 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++
2878
2879 for (ii = 0; ii < inum; ii++)
2880 {
2881 i = ilist[ii] ; itype = map[type[i]] ;
2882
2883 xtmp = x[i][0];
2884 ytmp = x[i][1];
2885 ztmp = x[i][2];
2886
2887
2888 jlist = firstneigh[i];
2889 jnum = numneigh[i];
2890 for (jj = 0; jj < jnum; jj++ )
2891 {
2892 j = jlist[jj] ;
2893 j &= NEIGHMASK;
2894 jtype = map[type[j]];
2895 if (jtype == itype) continue;
2896 m = intype[itype][jtype];
2897
2898 delr[0] = x[j][0] - xtmp;
2899 delr[1] = x[j][1] - ytmp;
2900 delr[2] = x[j][2] - ztmp;
2901 rsq = dot3(delr,delr);
2902
2903 if (sqrt(rsq) <= intparams[m].dc2) {
2904 flag_QEq[i] = 1; flag_QEq[j] = 1;
2905 }
2906 }
2907 if (flag_QEq[i] == 1) {
2908 QEq = 1;
2909 }
2910 }
2911
2912 // ::::::::::::::: Testing the presence of oxide :::::::::::
2913 m = 0;
2914 MPI_Allgather(&QEq,1,MPI_INT,QEqall,1,MPI_INT,world);
2915 for (iproc = 0; iproc < nproc; iproc++) {
2916 if (QEqall[iproc] == 1) m = 1;
2917 }
2918 if (m == 0) {
2919 memory->destroy(tabtemp);
2920 memory->destroy(Alltabtemp);
2921 memory->destroy(gptmp);
2922 memory->destroy(Allgptmp);
2923 memory->destroy(flag_gp);
2924 memory->destroy(tab_gp);
2925 memory->destroy(nelt);
2926
2927 memory->destroy(team_elt);
2928 memory->destroy(team_QEq);
2929 memory->destroy(QEqall);
2930 return;
2931 }
2932 // :::::::::::::::::::::::::::::::::::::::::::::::::::::::
2933
2934
2935 for (m = 0; m < nproc; m++) {
2936 for (i = 0; i < nall; i++) { flag_gp[m][i] = 0; }
2937 }
2938
2939 // OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
2940 // It includes oxygens entering the QEq scheme O
2941 // OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
2942
2943
2944 ngp = igp = 0; nelt[ngp] = 0;
2945
2946 // On prend un oxygène
2947 // printf ("[me %d] On prend un oxygene\n",me);
2948
2949 for (ii = 0; ii < inum; ii++) {
2950 i = ilist[ii] ; itype = map[type[i]];
2951 if (itype != 0 || flag_QEq[i] == 0) continue;
2952
2953 m = 0;
2954
2955 if (ngp != 0 && flag_gp[me][i] == ngp) continue;
2956
2957
2958 // Grouping Initialisation
2959 // -----------------------------
2960 if (flag_gp[me][i] == 0) {
2961 ngp += 1; nelt[ngp] = 0;
2962 tab_gp[ngp][nelt[ngp]] = i;
2963 flag_gp[me][i] = ngp;
2964 nelt[ngp] += 1;
2965 }
2966 // -------------------------------
2967
2968
2969 // Loop on the groups
2970 // ----------------------
2971 for (kk = 0; kk < nelt[ngp]; kk++)
2972 {
2973 k = tab_gp[ngp][kk];
2974 ktype = map[type[k]];
2975 // printf ("[me %d] kk - gp %d elemt %d : atom %d(%d)\n",me,ngp,kk,k,ktype);
2976 if (k >= nlocal) continue;
2977
2978 xtmp = x[k][0];
2979 ytmp = x[k][1];
2980 ztmp = x[k][2];
2981
2982 // Loop on the oxygen's neighbor of the group
2983 // ---------------------------------------------
2984 jlist = firstneigh[k];
2985 jnum = numneigh[k];
2986 for (j = 0; j < nall; j++ )
2987 {
2988 jtype = map[type[j]];
2989 if (jtype == ktype) continue;
2990 m = intype[itype][jtype];
2991
2992 if (jtype == 0 && flag_QEq[j] == 0) continue;
2993
2994
2995 if (flag_gp[me][j] == ngp) continue;
2996
2997 delr[0] = x[j][0] - xtmp;
2998 delr[1] = x[j][1] - ytmp;
2999 delr[2] = x[j][2] - ztmp;
3000 rsq = dot3(delr,delr);
3001
3002 // -------------------------------------
3003 if (sqrt(rsq) <= cutmax) {
3004
3005 flag_QEq[j] = 1; //Entre dans le schema QEq
3006
3007 // :::::::::::::::::::: Meeting of two group in the same proc :::::::::::::::::::::
3008
3009 if (flag_gp[me][j] != 0 && flag_gp[me][j] != ngp && nelt[flag_gp[me][j]] != 0) {
3010 printf("[me %d] (atom %d) %d [elt %d] rencontre un nouveau groupe %d [elt %d] (atom %d)\n",
3011 me,k,ngp,nelt[ngp],flag_gp[me][j],nelt[flag_gp[me][j]],j);
3012
3013 // On met a jours les tableaux
3014 // -----------------------------
3015 igp = flag_gp[me][j];
3016 z = min(igp,ngp);
3017
3018 if (z == igp) { igp = z; }
3019 else if (z == ngp) {
3020 ngp = igp ; igp = z;
3021 flag_gp[me][j] = ngp;
3022 }
3023
3024 for (zz = 0; zz < nelt[ngp]; zz++) {
3025 z = tab_gp[ngp][zz];
3026 tab_gp[igp][nelt[igp]] = z;
3027 nelt[igp] += 1;
3028 flag_gp[me][z] = igp;
3029 tab_gp[ngp][zz] = 0;
3030 }
3031
3032 nelt[ngp] = 0;
3033 for (z = nlocal; z < nall; z++) {
3034 if (flag_gp[me][z] == ngp) flag_gp[me][z] = igp;
3035 }
3036
3037 m = 1; kk = 0;
3038 ngp = igp;
3039 break;
3040 }
3041 // ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
3042
3043 flag_gp[me][j] = ngp;
3044 if (j < nlocal)
3045 {
3046 tab_gp[ngp][nelt[ngp]] = j;
3047 nelt[ngp] += 1;
3048 }
3049 }
3050 } // for j
3051 } // for k
3052 } // for ii
3053
3054 // OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
3055 // Groups communication
3056 // OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
3057 for (i = 0; i < nproc; i++) {
3058 forward_int(flag_gp[i]); //reverse_int(flag_gp[i]);
3059 }
3060 // ---
3061
3062
3063 // =======================================================
3064 // Loop on the cation to make them joined in the oxygen's
3065 // group which it interacts
3066 // =======================================================
3067 igp = 0;
3068 for (ii = 0; ii < inum; ii++)
3069 {
3070 i = ilist[ii] ; itype = map[type[i]];
3071 if (itype == 0) continue;
3072
3073 xtmp = x[i][0];
3074 ytmp = x[i][1];
3075 ztmp = x[i][2];
3076
3077 jlist = firstneigh[i];
3078 jnum = numneigh[i];
3079 for (jj = 0; jj < jnum; jj++ )
3080 {
3081 j = jlist[jj] ;
3082 j &= NEIGHMASK;
3083 jtype = map[type[j]];
3084 if (jtype != 0) continue;
3085
3086 m = 0;
3087 for (iproc = 0; iproc < nproc; iproc++) {
3088 if (flag_gp[iproc][j] != 0) m = flag_gp[iproc][j];
3089 }
3090 if (m == 0) continue;
3091
3092 delr[0] = x[j][0] - xtmp;
3093 delr[1] = x[j][1] - ytmp;
3094 delr[2] = x[j][2] - ztmp;
3095 rsq = dot3(delr,delr);
3096
3097 // ----------------------------------------
3098 if (sqrt(rsq) <= cutmax) {
3099 // if (sqrt(rsq) <= intparams[m].dc2) {
3100 // ----------------------------------------
3101
3102 flag_QEq[i] = 1; igp = flag_gp[me][j];
3103
3104 if (flag_gp[me][i] == 0) flag_gp[me][i] = igp;
3105
3106 if (flag_gp[me][i] != igp && igp != 0) {
3107 printf ("[me %d] Cation i %d gp %d [nelt %d] rencontre j %d(%d)du groupe %d [nelt %d]\n",
3108 me,i,flag_gp[me][i],nelt[flag_gp[me][i]],j,jtype,igp,nelt[igp]);
3109
3110 igp = min(flag_gp[me][i],flag_gp[me][j]);
3111 if (igp == flag_gp[me][i]) { kgp = flag_gp[me][j]; }
3112 else { kgp = flag_gp[me][i]; }
3113
3114 for (k = 0; k < nelt[kgp]; k++) {
3115 z = tab_gp[kgp][k];
3116 tab_gp[igp][nelt[igp]] = z;
3117 nelt[igp] += 1;
3118 flag_gp[me][z] = igp;
3119 tab_gp[kgp][k] = 0;
3120 }
3121 nelt[kgp] = 0;
3122
3123 for (k = 0; k < nall; k++) {
3124 if (flag_gp[me][k] == kgp) flag_gp[me][k] = igp;
3125 }
3126
3127 }
3128 m = 0;
3129 for (k = 0; k < nelt[igp]; k++) {
3130 if (tab_gp[igp][k] == i) m = 1;
3131 }
3132
3133 if (i >= nlocal || m == 1 ) continue;
3134 // printf ("[me %d] igp %d - nelt %d atom %d\n",me,igp,nelt[igp],i);
3135 tab_gp[igp][nelt[igp]] = i;
3136 nelt[igp] += 1;
3137 break;
3138 }
3139
3140 } // voisin j
3141
3142 } // atom i
3143
3144 /* ==================================================
3145 Group Communication between proc for unification
3146 ================================================== */
3147 for (i = 0; i < nproc; i++) {
3148 forward_int(flag_gp[i]);// reverse_int(flag_gp[i]);
3149 }
3150
3151 // =============== End of COMM =================
3152
3153
3154 for (i = 0; i < nall; i++) {
3155
3156 m = 10*me + flag_gp[me][i];
3157 if (m == 10*me) continue; // Pas de groupe zero
3158 gptmp[m] = 1;
3159 for (k = 0; k < nproc; k++) {
3160
3161 if (k == me) continue;
3162 if (tabtemp[m][k] != 0) continue;
3163
3164 if (flag_gp[k][i] != 0) {
3165 tabtemp[m][k] = 10*k + flag_gp[k][i];
3166 }
3167 }
3168
3169 }
3170
3171 for (k = 0; k < 10*nproc; k++) {
3172 MPI_Allreduce(tabtemp[k],Alltabtemp[k],nproc,MPI_INT,MPI_SUM,world); }
3173 MPI_Allreduce(gptmp,Allgptmp,10*nproc,MPI_INT,MPI_SUM,world);
3174
3175 nteam = 0; iproc = 0;
3176 for (igp = 0; igp < 10*nproc; igp++) {
3177 if (Allgptmp[igp] == 0) continue;
3178 iproc = int(static_cast<double>(igp)/10.0);
3179 ngp = igp - 10*iproc;
3180 if (nteam == 0) {
3181
3182 nteam += 1;
3183 team_elt[nteam][iproc] = 0;
3184 team_QEq[nteam][iproc][team_elt[nteam][iproc]] = ngp;
3185 team_elt[nteam][iproc] += 1;
3186 } else {
3187 m = 0;
3188 for (i = 1; i < nteam+1; i++) {
3189 for (k = 0; k < team_elt[i][iproc]; k++) {
3190 if (ngp == team_QEq[i][iproc][k]) m = 1;
3191 } }
3192 if (m == 1) continue;
3193 // create a new team!!
3194 // ---------------------------
3195 if (m == 0) {
3196 nteam += 1;
3197 team_elt[nteam][iproc] = 0;
3198 team_QEq[nteam][iproc][team_elt[nteam][iproc]] = ngp;
3199 team_elt[nteam][iproc] += 1;
3200 }
3201 }
3202 // -------
3203 // On a mis une graine dans le groupe et nous allons
3204 // remplir se groupe en questionnant le "tabtemp[igp][iproc]=jgp"
3205 //
3206 for (kk = 0; kk < nproc; kk++) {
3207 for (k = 0; k < team_elt[nteam][kk]; k++) {
3208
3209 // On prend le gp le k ieme element de la team nteam sur le proc iproc
3210 // ngp = 0;
3211 ngp = team_QEq[nteam][kk][k];
3212 kgp = 10*kk + ngp;
3213
3214 // On regarde sur les autre proc si ce gp ne pointe pas vers un autre.
3215 for (i = 0; i < nproc; i++) {
3216 if (i == kk) continue;
3217 if (Alltabtemp[kgp][i] == 0) continue;
3218
3219 if (Alltabtemp[kgp][i] != 0) ngp = Alltabtemp[kgp][i];
3220 ngp = ngp - 10*i;
3221
3222 // Est ce que ce groupe est nouveau?
3223 m = 0;
3224 for (j = 0; j < team_elt[nteam][i]; j++) {
3225 if (team_QEq[nteam][i][j] == ngp) m = 1;
3226 }
3227
3228 if (m == 0) {
3229 iproc = i; k = 0;
3230 team_QEq[nteam][i][team_elt[nteam][i]] = ngp ;
3231 team_elt[nteam][i] += 1;
3232 }
3233 } // regard sur les autre proc
3234
3235 } // On rempli de proche en proche
3236 } // boucle kk sur les proc
3237 }
3238
3239 // Finalement on met le numero de la team en indice du flag_QEq, c mieu!
3240 // ---------------------------------------------------------------------
3241
3242 for (ii = 0; ii < inum; ii++)
3243 {
3244 i = ilist[ii]; m = 0; itype = map[type[i]];
3245 if (flag_QEq[i] == 0) continue;
3246
3247 gp = flag_gp[me][i];
3248 for (j = 1; j <= nteam; j++) {
3249 for (k = 0; k < team_elt[j][me]; k++) {
3250 if (gp == team_QEq[j][me][k]) {
3251 flag_QEq[i] = j; m = 1;
3252 break;
3253 }
3254 }
3255 if (m == 1) break;
3256 }
3257 }
3258
3259 memory->destroy(tabtemp);
3260 memory->destroy(Alltabtemp);
3261 memory->destroy(gptmp);
3262 memory->destroy(Allgptmp);
3263 memory->destroy(flag_gp);
3264 memory->destroy(tab_gp);
3265 memory->destroy(nelt);
3266
3267 memory->destroy(team_elt);
3268 memory->destroy(team_QEq);
3269 memory->destroy(QEqall);
3270 }
3271
3272 /* ---------------------------------------------------------------------- */
3273
Init_charge(int *,int *,int *)3274 void PairSMTBQ::Init_charge(int * /*nQEq*/, int * /*nQEqa*/, int * /*nQEqc*/)
3275 {
3276 int ii,i,gp,itype;
3277 int *ilist;
3278 std::vector<int> test(cluster),init(cluster);
3279 double bound,tot,totll;
3280
3281 int inum = list->inum;
3282 int *type = atom->type;
3283 double *q = atom->q;
3284 ilist = list->ilist;
3285
3286 if (nteam == 0) return;
3287
3288 if (me == 0) printf (" ======== Init_charge ======== \n");
3289 for (gp = 0; gp < cluster; gp++) {
3290 test[gp] = 0; init[gp] = 0;
3291 }
3292
3293
3294 // On fait un test sur les charges pour voir sont
3295 // elles sont dans le domaine delimiter par DeltaQ
3296 // -------------------------------------------------
3297 for (ii = 0; ii < inum; ii++) {
3298 i = ilist[ii]; itype = map[type[i]];
3299 gp = flag_QEq[i];
3300
3301 if (gp != 0 && itype == 0) {
3302 bound = fabs(2.0*ncov/params[itype].sto - fabs(params[itype].qform)) ;
3303 if (bound == fabs(params[itype].qform)) continue;
3304 if (fabs(q[i]) < bound) test[gp] = 1;
3305 }
3306 }
3307
3308 // TODO
3309 MPI_Allreduce(test.data(),init.data(),cluster,MPI_INT,MPI_SUM,world);
3310
3311 // On fait que sur les atomes hybrides!!!
3312 // ----------------------------------------
3313
3314 tot = totll = 0.0;
3315 for (ii = 0; ii < inum; ii++) {
3316 i = ilist[ii]; itype = map[type[i]];
3317 gp = flag_QEq[i];
3318
3319 if (gp != 0 && init[gp] != 0) {
3320 if (itype == 0) q[i] = -1.96;
3321 if (itype != 0) q[i] = 1.96 * static_cast<double>(nQEqaall[gp]) / static_cast<double>(nQEqcall[gp]);
3322 }
3323 tot += q[i];
3324 }
3325 MPI_Allreduce(&tot,&totll,1,MPI_INT,MPI_SUM,world);
3326 if (me == 0) printf (" === Fin de init_charge qtot %20.15f ====\n",totll);
3327
3328 }
3329 /* ----------------------------------------------------------------------
3330 * COMMUNICATION
3331 * ---------------------------------------------------------------------- */
3332
pack_forward_comm(int n,int * list,double * buf,int,int *)3333 int PairSMTBQ::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
3334 {
3335 int i,j,m;
3336
3337 m = 0;
3338 for (i = 0; i < n; i ++) {
3339 j = list[i];
3340 buf[m++] = tab_comm[j];
3341 // if (j < 3) printf ("[%d] %d pfc %d %d buf_send = %f \n",me,n,i,m-1,buf[m-1]);
3342 }
3343 return m;
3344 }
3345
3346 /* ---------------------------------------------------------------------- */
3347
unpack_forward_comm(int n,int first,double * buf)3348 void PairSMTBQ::unpack_forward_comm(int n, int first, double *buf)
3349 {
3350 int i,m,last;
3351
3352 m = 0;
3353 last = first + n ;
3354 for (i = first; i < last; i++) {
3355 tab_comm[i] = buf[m++];
3356 // if (i<first+3) printf ("[%d] ufc %d %d buf_recv = %f \n",me,i,m-1,buf[m-1]);
3357 }
3358 }
3359
3360 /* ---------------------------------------------------------------------- */
3361
pack_reverse_comm(int n,int first,double * buf)3362 int PairSMTBQ::pack_reverse_comm(int n, int first, double *buf)
3363 {
3364 int i,m,last;
3365
3366 m = 0;
3367 last = first + n;
3368 for (i = first; i < last; i++) {
3369 buf[m++] = tab_comm[i];
3370 // if (i<first+3) printf ("[%d] prc %d %d buf_send = %f \n",me,i,m-1,buf[m-1]);
3371 }
3372 return m;
3373 }
3374
3375 /* ---------------------------------------------------------------------- */
3376
unpack_reverse_comm(int n,int * list,double * buf)3377 void PairSMTBQ::unpack_reverse_comm(int n, int *list, double *buf)
3378 {
3379 int i,j,m;
3380
3381 m = 0;
3382 for (i = 0; i < n; i++) {
3383 j = list[i];
3384 // tab_comm[j] += buf[m++];
3385 tab_comm[j] = buf[m++];
3386 // if (j<3) printf ("[%d] %d urc %d %d buf_recv = %f \n",me,n,i,m-1,buf[m-1]);
3387 }
3388 }
3389
3390 /* ---------------------------------------------------------------------- */
3391
forward(double * tab)3392 void PairSMTBQ::forward(double *tab)
3393 {
3394 int i;
3395 int nlocal = atom->nlocal;
3396 int nghost = atom->nghost;
3397
3398 for (i=0; i<nlocal+nghost; i++) tab_comm[i] = tab[i];
3399
3400 comm->forward_comm_pair(this);
3401
3402 for (i=0; i<nlocal+nghost; i++) tab[i] = tab_comm[i];
3403 }
3404
3405 /* ---------------------------------------------------------------------- */
3406
reverse(double * tab)3407 void PairSMTBQ::reverse(double *tab)
3408 {
3409 int i;
3410 int nlocal = atom->nlocal;
3411 int nghost = atom->nghost;
3412
3413 for (i=0; i<nlocal+nghost; i++) tab_comm[i] = tab[i];
3414
3415 comm->reverse_comm_pair(this);
3416
3417 for (i=0; i<nlocal+nghost; i++) tab[i] = tab_comm[i];
3418 }
3419
3420 /* ---------------------------------------------------------------------- */
3421
forward_int(int * tab)3422 void PairSMTBQ::forward_int(int *tab)
3423 {
3424 int i;
3425 int nlocal = atom->nlocal;
3426 int nghost = atom->nghost;
3427
3428 for (i=0; i<nlocal+nghost; i++) { tab_comm[i] = static_cast<double>(tab[i]);}
3429
3430 comm->forward_comm_pair(this);
3431
3432 for (i=0; i<nlocal+nghost; i++) {
3433 if (fabs(tab_comm[i]) > 0.1) tab[i] = int(tab_comm[i]) ; }
3434 }
3435
3436 /* ---------------------------------------------------------------------- */
3437
reverse_int(int * tab)3438 void PairSMTBQ::reverse_int(int *tab)
3439 {
3440 int i;
3441 int nlocal = atom->nlocal;
3442 int nghost = atom->nghost;
3443
3444 for (i=0; i<nlocal+nghost; i++) { tab_comm[i] = static_cast<double>(tab[i]);}
3445
3446 comm->reverse_comm_pair(this);
3447
3448 for (i=0; i<nlocal+nghost; i++) {
3449 if (fabs(tab_comm[i]) > 0.1) tab[i] = int(tab_comm[i]); }
3450 }
3451
3452 /* ---------------------------------------------------------------------- */
3453 /* ---------------------------------------------------------------------- */
3454
memory_usage()3455 double PairSMTBQ::memory_usage()
3456 {
3457 double bytes = (double)maxeatom * sizeof(double);
3458 bytes += (double)maxvatom*6 * sizeof(double);
3459 bytes += (double)nmax * sizeof(int);
3460 bytes += (double)MAXNEIGH * nmax * sizeof(int);
3461 return bytes;
3462 }
3463
3464 /* ---------------------------------------------------------------------- */
3465
add_pages(int howmany)3466 void PairSMTBQ::add_pages(int howmany)
3467 {
3468 int toppage = maxpage;
3469 maxpage += howmany*PGDELTA;
3470
3471 pages = (int **)
3472 memory->srealloc(pages,maxpage*sizeof(int *),"pair:pages");
3473 for (int i = toppage; i < maxpage; i++)
3474 memory->create(pages[i],pgsize,"pair:pages[i]");
3475 }
3476
3477 /* ---------------------------------------------------------------------- */
3478
3479
Tokenize(char * s,char *** tok)3480 int PairSMTBQ::Tokenize( char* s, char*** tok )
3481 {
3482 char test[MAXLINE];
3483 const char *sep = "' ";
3484 char *mot;
3485 int count=0;
3486 mot = nullptr;
3487
3488
3489 strncpy( test, s, MAXLINE-1 );
3490
3491 for (mot = strtok(test, sep); mot; mot = strtok(nullptr, sep)) {
3492 strncpy( (*tok)[count], mot, MAXLINE );
3493 count++;
3494 }
3495
3496 return count;
3497 }
3498
3499
3500
CheckEnergyVSForce()3501 void PairSMTBQ::CheckEnergyVSForce()
3502 {
3503 double drL,iq,jq,rsq,evdwlCoul,fpairCoul,eflag=0,ErepR,frepR,fpair,evdwl;
3504 int i,j,iiiMax,iii,iCoord;
3505 int itype,jtype,l,m;
3506 double r,t1,t2,sds,xi,engSurf,fforceSurf;
3507 double eng,fforce,engBB,fforceBB;
3508
3509 double za,zb,gam,dgam,dza,dzb,
3510 d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2,na,nb;
3511 int *type = atom->type;
3512 const char *NameFile;
3513
3514
3515 i=0;
3516 j=1;
3517 map[type[i]]=0; //ox
3518 itype=map[type[i]];
3519 iq=-1;
3520
3521
3522 map[type[j]]=0; //ox
3523 jtype=map[type[j]];
3524 jq=-1;
3525 coord[i]=coordOxBulk;
3526 coord[j]=coordOxBulk;
3527 m = intype[itype][jtype];
3528
3529
3530 na = params[itype].ne ;
3531 nb = params[jtype].ne ;
3532 za = params[itype].dzeta ;
3533 zb = params[jtype].dzeta ;
3534
3535
3536 // Ouverture du fichier
3537
3538 for (iCoord=1;iCoord<5; iCoord++)
3539 {
3540
3541 if (iCoord==1)
3542 {coord[i]=2.2;
3543 coord[j]=2.1;
3544 NameFile=(const char *)"energyandForceOxOxUnderCoord.txt";
3545 }
3546 if (iCoord==2)
3547 {coord[i]=coordOxBulk;
3548 coord[j]=coordOxBulk;
3549 NameFile=(const char *)"energyandForceOxOxCoordBulk.txt";
3550 }
3551 if (iCoord==3)
3552 {coord[i]=3.8;
3553 coord[j]=4;
3554 NameFile=(const char *)"energyandForceOxOxOverCoord.txt";
3555 }
3556 if (iCoord==4)
3557 {coord[i]=2.2;
3558 coord[j]=3.5;
3559 NameFile=(const char *)"energyandForceOxOxUnderOverCoord.txt";
3560 }
3561
3562
3563 ofstream fichierOxOx(NameFile, ios::out | ios::trunc) ;
3564
3565 drL=0.0001;
3566 iiiMax=int((cutmax-1.2)/drL);
3567 for (iii=1; iii< iiiMax ; iii++) {
3568 r=1.2+drL*iii;
3569 rsq=r*r;
3570 evdwlCoul = 0.0 ; fpairCoul = 0.0;
3571 potqeq(i,j,iq,jq,rsq,fpairCoul,eflag,evdwlCoul);
3572 fpairCoul=fpairCoul*r;
3573
3574 rep_OO (&intparams[m],rsq,fpair,eflag,evdwl);
3575 ErepR = evdwl;
3576 frepR= fpair*r;
3577
3578 gam = dgam = dza = dzb = d2zaa = d2zab =
3579 d2zbb = d2zra = d2zrb = d2gamr2 = 0.0 ;
3580
3581
3582 // gammas_(na,nb,za,zb,r,gam,dgam,dza,dzb,
3583 // d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
3584 gammas(na,nb,za,zb,r,gam,dgam,dza,dzb,
3585 d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
3586
3587
3588 sds = rsq/ds ; l = int(sds) ;
3589 xi = sds - static_cast<double>(l) ;
3590
3591
3592 t1 = fafb[l][m] + (fafb[l+1][m] - fafb[l][m])*xi;
3593 t2 = fafb[l+1][m] + (fafb[l+2][m] - fafb[l+1][m])*(xi-1.0);
3594 eng = iq*jq*(t1 + (t2 - t1)*xi/2.0);
3595
3596 t1 = dfafb[l][m] + (dfafb[l+1][m] - dfafb[l][m])*xi;
3597 t2 = dfafb[l+1][m] + (dfafb[l+2][m] - dfafb[l+1][m])*(xi-1);
3598 fforce = - iq*jq*(t1 + (t2 - t1)*xi/2.0)*r ;
3599
3600
3601 t1 = fafbOxOxSurf[l] + (fafbOxOxSurf[l+1] - fafbOxOxSurf[l])*xi;
3602 t2 = fafbOxOxSurf[l+1] + (fafbOxOxSurf[l+2] - fafbOxOxSurf[l+1])*(xi-1.0);
3603 engSurf = iq*jq*(t1 + (t2 - t1)*xi/2.0);
3604
3605 t1 = fafbOxOxBB[l] + (fafbOxOxBB[l+1] - fafbOxOxBB[l])*xi;
3606 t2 = fafbOxOxBB[l+1] + (fafbOxOxBB[l+2] - fafbOxOxBB[l+1])*(xi-1.0);
3607 engBB = iq*jq*(t1 + (t2 - t1)*xi/2.0);
3608
3609 t1 = dfafbOxOxSurf[l] + (dfafbOxOxSurf[l+1] - dfafbOxOxSurf[l])*xi;
3610 t2 = dfafbOxOxSurf[l+1] + (dfafbOxOxSurf[l+2] - dfafbOxOxSurf[l+1])*(xi-1);
3611 fforceSurf = - iq*jq*(t1 + (t2 - t1)*xi/2.0)*r ;
3612
3613 t1 = dfafbOxOxBB[l] + (dfafbOxOxBB[l+1] - dfafbOxOxBB[l])*xi;
3614 t2 = dfafbOxOxBB[l+1] + (dfafbOxOxBB[l+2] - dfafbOxOxBB[l+1])*(xi-1);
3615 fforceBB = - iq*jq*(t1 + (t2 - t1)*xi/2.0)*r ;
3616
3617 if (fichierOxOx) { fichierOxOx<< setprecision (9) <<r <<" "<<evdwlCoul <<" " <<fpairCoul <<" "<<eng <<" " <<fforce <<" "<<engSurf <<" " <<fforceSurf <<" "<<engBB <<" " <<fforceBB <<" "<<ErepR<<" "<<frepR<<" "<<gam<<" "<<dgam<<endl ;}
3618
3619 }
3620
3621
3622 if (fichierOxOx) fichierOxOx.close() ;
3623 }
3624
3625
3626
3627 map[type[j]]=1; //met
3628 jtype=map[type[j]];
3629 jq=1;
3630 coord[i]=coordOxBulk;
3631 coord[j]=6;
3632 m = intype[itype][jtype];
3633
3634
3635 na = params[itype].ne ;
3636 nb = params[jtype].ne ;
3637 za = params[itype].dzeta ;
3638 zb = params[jtype].dzeta ;
3639
3640
3641 // Ouverture du fichier
3642
3643 for (iCoord=1;iCoord<4; iCoord++)
3644 {
3645
3646 if (iCoord==1)
3647 {coord[i]=2.2;
3648 coord[j]=2.1;
3649 NameFile="energyandForceOxTiUnderCoord.txt";
3650 }
3651 if (iCoord==2)
3652 {coord[i]=coordOxBulk;
3653 coord[j]=coordOxBulk;
3654 NameFile="energyandForceOxTiCoordBulk.txt";
3655 }
3656 if (iCoord==3)
3657 {coord[i]=3.8;
3658 coord[j]=4;
3659 NameFile="energyandForceOxTiOverCoord.txt";
3660 }
3661
3662
3663 ofstream fichierOxTi(NameFile, ios::out | ios::trunc) ;
3664
3665 drL=0.0001;
3666 iiiMax=int((cutmax-1.2)/drL);
3667 for (iii=1; iii< iiiMax ; iii++) {
3668 r=1.2+drL*iii;
3669 rsq=r*r;
3670 evdwlCoul = 0.0 ; fpairCoul = 0.0;
3671 potqeq(i,j,iq,jq,rsq,fpairCoul,eflag,evdwlCoul);
3672 fpairCoul=fpairCoul*r;
3673
3674 rep_OO (&intparams[m],rsq,fpair,eflag,evdwl);
3675 ErepR = evdwl;
3676 frepR= fpair*r;
3677
3678 gam = dgam = dza = dzb = d2zaa = d2zab =
3679 d2zbb = d2zra = d2zrb = d2gamr2 = 0.0 ;
3680
3681
3682 // gammas_(na,nb,za,zb,r,gam,dgam,dza,dzb,
3683 // d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
3684 gammas(na,nb,za,zb,r,gam,dgam,dza,dzb,
3685 d2zaa,d2zab,d2zbb,d2zra,d2zrb,d2gamr2) ;
3686
3687
3688 if (fichierOxTi) { fichierOxTi<< setprecision (9) <<r <<" "<<evdwlCoul <<" " <<fpairCoul <<" "<<ErepR<<" "<<frepR<<" "<<gam<<" "<<dgam<<endl ;}
3689
3690 }
3691
3692
3693 if (fichierOxTi) fichierOxTi.close() ;
3694 }
3695
3696 exit(0);
3697 }
3698
3699 /* :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
3700 GAMMAS FUNCTION (GALE)
3701 ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: */
3702
gammas(double & na,double & nb,double & za,double & zb,double & r,double & gam,double & dgam,double & dza,double & dzb,double & d2zaa,double & d2zab,double & d2zbb,double & d2zra,double & d2zrb,double & d2gamr2)3703 void PairSMTBQ::gammas(double &na, double &nb, double &za, double &zb, double &r, double &gam,
3704 double &dgam, double &dza, double &dzb, double &d2zaa, double &d2zab, double &d2zbb,
3705 double &d2zra, double &d2zrb, double &d2gamr2)
3706 {
3707 /* ---------------------------------------------------------------
3708 Subroutine calculates the integral over two s orbtials
3709 for Slater functions
3710
3711 On input :
3712
3713 na = principle quantum number of A
3714 nb = principle quantum number of B
3715 za = orbital exponent of A
3716 zb = orbital exponent of B
3717 r = distance between A and B
3718
3719 On exit :
3720
3721 gam = integral
3722 dgam = first derivative of gam with respect to r
3723 dza = first derivative of gam with respect to za
3724 dzb = first derivative of gam with respect to zb
3725 d2zaa = second derivative of gam with respect to za/za
3726 d2zab = second derivative of gam with respect to za/zb
3727 d2zbb = second derivative of gam with respect to zb/zb
3728 d2zra = second derivative of gam with respect to r/za
3729 d2zrb = second derivative of gam with respect to r/zb
3730 d2gamr2 = second derivative of gam with respect to r
3731
3732 Julian Gale, Imperial College, December 1997
3733 ---------------------------------------------------------------- */
3734
3735 int i;
3736 double z2ra,z2rb,na2,nb2,halfr,drtrm,rtrm,ss,deriv,
3737 dzeta1,dzeta2,d2zeta11,d2zeta12,d2zeta22,d2zeta1r,
3738 d2zeta2r,deriv2,trm1,trm2,trm3,d2rtrm,ztrm,ctrm,
3739 rfct1,rgam1,rdza1,rgam2;
3740
3741
3742
3743 gam=0.0;
3744 dgam=0.0;
3745 dza=0.0;
3746 dzb=0.0;
3747 d2zaa=0.0;
3748 d2zab=0.0;
3749 d2zbb=0.0;
3750 d2zra=0.0;
3751 d2zrb=0.0;
3752 d2gamr2=0.0;
3753
3754 // This routine only handles two centre integrals
3755
3756
3757 if (r < 1.0e-10) return;
3758
3759
3760 // Create local variables
3761
3762 z2ra=2.0*za*r;
3763 z2rb=2.0*zb*r;
3764 na2=2*na;
3765 nb2=2*nb;
3766 halfr=0.5*r;
3767 d2rtrm=powint(halfr,na2-2);
3768 drtrm=d2rtrm*halfr;
3769 rtrm=drtrm*halfr;
3770
3771 // First term
3772
3773
3774 css(ss,na2-1,0,z2ra,0.0,r,deriv,dzeta1,dzeta2,
3775 d2zeta11,d2zeta12,d2zeta22,d2zeta1r,d2zeta2r,
3776 deriv2) ;
3777
3778 gam=rtrm*ss;
3779 dgam=rtrm*deriv+static_cast<double>(na)*drtrm*ss;
3780 dza=rtrm*dzeta1;
3781 d2zaa=rtrm*d2zeta11;
3782 d2zra=rtrm*d2zeta1r+static_cast<double>(na)*drtrm*dzeta1;
3783 d2gamr2=d2gamr2+0.5*static_cast<double>(na*(na2-1))*d2rtrm*ss + 2.0*static_cast<double>(na)*drtrm*deriv+rtrm*deriv2;
3784
3785 // Sum over 2*nb
3786
3787 rtrm=drtrm;
3788 drtrm=d2rtrm;
3789 ztrm=0.5/(zb*static_cast<double>(nb2));
3790
3791 for (i = nb2; i >= 1; i--) {
3792 rtrm=rtrm*halfr;
3793 drtrm=drtrm*halfr;
3794 ztrm=ztrm*2.0*zb;
3795 ctrm=ztrm/factorial(static_cast<int>(nb2-i));
3796
3797 css(ss,na2-1,nb2-i,z2ra,z2rb,r,deriv,dzeta1,dzeta2,
3798 d2zeta11,d2zeta12,d2zeta22,d2zeta1r,d2zeta2r,deriv2);
3799
3800 trm1=static_cast<double>(i)*ctrm;
3801 trm2=trm1*rtrm;
3802 gam=gam-trm2*ss;
3803 trm3=trm1*static_cast<double>(na2+nb2-i)*drtrm;
3804 dgam=dgam-trm2*deriv-0.5*trm3*ss;
3805 d2gamr2=d2gamr2-trm2*deriv2-trm3*deriv-0.5*trm3*static_cast<double>(na2+nb2-i-1)*ss/r;
3806 dza=dza-trm2*dzeta1;
3807 dzb=dzb-(trm2/zb)*((static_cast<double>(nb2-i))*ss+zb*dzeta2);
3808 d2zaa=d2zaa-trm2*d2zeta11;
3809 d2zab=d2zab-(trm2/zb)*((static_cast<double>(nb2-i))*dzeta1+zb*d2zeta12);
3810 d2zbb=d2zbb-(trm2/zb)*(2.0*(static_cast<double>(nb2-i))*dzeta2+zb*d2zeta22 +
3811 (static_cast<double>((nb2-i-1)*(nb2-i))*ss/zb));
3812 d2zra=d2zra-trm2*d2zeta1r-0.5*trm3*dzeta1;
3813 d2zrb=d2zrb-(trm2/zb)*((static_cast<double>(nb2-i))*deriv+zb*d2zeta2r) -
3814 0.5*(trm3/zb)*((static_cast<double>(nb2-i))*ss+zb*dzeta2);
3815 }
3816
3817 // Multiply by coefficients
3818
3819 trm3=powint(2.0*za,na2+1)/factorial(static_cast<int>(na2));
3820 gam=gam*trm3;
3821 dgam=dgam*trm3;
3822 rfct1=((static_cast<double>(na2+1))/za);
3823 rgam1=rfct1*gam;
3824 dza=dza*trm3;
3825 rdza1=2.0*rfct1*dza;
3826 dza=dza+rgam1;
3827 dzb=dzb*trm3;
3828 rgam2=rgam1*static_cast<double>(na2)/za;
3829 d2zaa=d2zaa*trm3+rgam2+rdza1;
3830 d2zab=d2zab*trm3+rfct1*dzb;
3831 d2zbb=d2zbb*trm3;
3832 d2zra=d2zra*trm3+rfct1*dgam;
3833 d2zrb=d2zrb*trm3;
3834 d2gamr2=d2gamr2*trm3;
3835 return;
3836
3837 }
3838 /* --------------------------------------------------------------------------------
3839 Css
3840 -------------------------------------------------------------------------------- */
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)3841 void PairSMTBQ::css(double &s, double nn1, double nn2, double alpha, double beta, double r,
3842 double &deriv, double &dzeta1, double &dzeta2, double &d2zeta11, double &d2zeta12,
3843 double &d2zeta22, double &d2zeta1r, double &d2zeta2r, double &deriv2)
3844 {
3845 // implicit real (a-h,o-z)
3846 // common /fctrl/ fct(30) // A RAJOUTER DANS Pair_SMTBQ.h
3847 /* ------------------------------------------------------------------
3848 Modified integral calculation routine for Slater orbitals
3849 including derivatives. This version is for S orbitals only.
3850
3851 dzeta1 and dzeta2 are the first derivatives with respect to zetas
3852 and d2zeta11/d2zeta12/d2zeta22 are the second.
3853 d2zeta1r and d2zeta2r are the mixed zeta/r second derivatives
3854 deriv2 is the second derivative with respect to r
3855
3856 Julian Gale, Imperial College, December 1997
3857 ------------------------------------------------------------------- */
3858 int i,i1,nni1; //ulim
3859 double ulim,n1,n2,p,pt,x,k,dpdz1,dpdz2,dptdz1,dptdz2,dpdr,dptdr,d2pdz1r,
3860 d2pdz2r,d2ptdz1r,d2ptdz2r,zeta1,zeta2,sumzeta,difzeta,coff;
3861
3862 double da1[30],da2[30],db1[30],db2[30];
3863 double d2a11[30],d2a12[30],d2a22[30],dar[30];
3864 double d2b11[30],d2b12[30],d2b22[30],dbr[30];
3865 double d2a1r[30],d2a2r[30],d2b1r[30],d2b2r[30];
3866 double d2ar2[30],d2br2[30];
3867
3868 double *a,*b;
3869 memory->create(a,31,"pair:a");
3870 memory->create(b,31,"pair:a");
3871
3872
3873 // Set up factorials - stored as factorial(n) in location(n+1)
3874
3875 for (i = 1; i <= 30; i++) {
3876 fct[i]=factorial(i-1);
3877 }
3878 dzeta1=0.0;
3879 dzeta2=0.0;
3880 d2zeta11=0.0;
3881 d2zeta12=0.0;
3882 d2zeta22=0.0;
3883 d2zeta1r=0.0;
3884 d2zeta2r=0.0;
3885 deriv=0.0;
3886 deriv2=0.0;
3887 n1=nn1;
3888 n2=nn2;
3889 p =(alpha + beta)*0.5;
3890 pt=(alpha - beta)*0.5;
3891 x = 0.0;
3892 zeta1=alpha/r;
3893 zeta2=beta/r;
3894 sumzeta=zeta1+zeta2;
3895 difzeta=zeta1-zeta2;
3896
3897 // Partial derivative terms for zeta derivatives
3898
3899 dpdz1=r;
3900 dpdz2=r;
3901 dptdz1=r;
3902 dptdz2=-r;
3903 dpdr=0.5*sumzeta;
3904 dptdr=0.5*difzeta;
3905 d2pdz1r=1.0;
3906 d2pdz2r=1.0;
3907 d2ptdz1r=1.0;
3908 d2ptdz2r=-1.0;
3909
3910 // Reverse quantum numbers if necessary -
3911 // also change the sign of difzeta to match
3912 // change in sign of pt
3913
3914 if (n2 < n1) {
3915 k = n1;
3916 n1= n2;
3917 n2= k;
3918 pt=-pt;
3919 difzeta=-difzeta;
3920 dptdr=-dptdr;
3921 dptdz1=-dptdz1;
3922 dptdz2=-dptdz2;
3923 d2ptdz1r=-d2ptdz1r;
3924 d2ptdz2r=-d2ptdz2r;
3925 }
3926
3927 // Trap for enormously long distances which would cause
3928 // caintgs or cbintgs to crash with an overflow
3929
3930 if (p > 86.0 || pt > 86.0) {
3931 s=0.0;
3932 return;
3933 }
3934 //***************************
3935 // Find a and b integrals *
3936 //***************************
3937 caintgs(p,n1+n2+3,a);
3938 cbintgs(pt,n1+n2+3,b);
3939
3940 // Convert derivatives with respect to p and pt
3941 // into derivatives with respect to zeta1 and
3942 // zeta2
3943
3944 ulim=n1+n2+1;
3945 for (i = 1; i <= int(ulim); i++) {
3946 da1[i]=-a[i+1]*dpdz1;
3947 da2[i]=-a[i+1]*dpdz2;
3948 db1[i]=-b[i+1]*dptdz1;
3949 db2[i]=-b[i+1]*dptdz2;
3950 d2a11[i]=a[i+2]*dpdz1*dpdz1;
3951 d2a12[i]=a[i+2]*dpdz1*dpdz2;
3952 d2a22[i]=a[i+2]*dpdz2*dpdz2;
3953 d2b11[i]=b[i+2]*dptdz1*dptdz1;
3954 d2b12[i]=b[i+2]*dptdz1*dptdz2;
3955 d2b22[i]=b[i+2]*dptdz2*dptdz2;
3956 dar[i]=-a[i+1]*dpdr;
3957 dbr[i]=-b[i+1]*dptdr;
3958 d2a1r[i]=a[i+2]*dpdz1*dpdr-a[i+1]*d2pdz1r;
3959 d2a2r[i]=a[i+2]*dpdz2*dpdr-a[i+1]*d2pdz2r;
3960 d2b1r[i]=b[i+2]*dptdz1*dptdr-b[i+1]*d2ptdz1r;
3961 d2b2r[i]=b[i+2]*dptdz2*dptdr-b[i+1]*d2ptdz2r;
3962 d2ar2[i]=a[i+2]*dpdr*dpdr;
3963 d2br2[i]=b[i+2]*dptdr*dptdr;
3964 }
3965
3966 // Begin section used for overlap integrals involving s functions
3967
3968 for (i1 = 1; i1 <= int(ulim); i1++) {
3969 nni1=n1+n2-i1+2;
3970 coff=coeffs(n1,n2,i1-1);
3971 deriv=deriv+coff*(dar[i1]*b[nni1]+a[i1]*dbr[nni1]);
3972 x=x+coff*a[i1]*b[nni1];
3973 dzeta1=dzeta1+coff*(da1[i1]*b[nni1]+a[i1]*db1[nni1]);
3974 dzeta2=dzeta2+coff*(da2[i1]*b[nni1]+a[i1]*db2[nni1]);
3975 d2zeta11=d2zeta11+coff*(d2a11[i1]*b[nni1]+a[i1]*d2b11[nni1]+
3976 2.0*da1[i1]*db1[nni1]);
3977 d2zeta12=d2zeta12+coff*(d2a12[i1]*b[nni1]+a[i1]*d2b12[nni1]+
3978 da1[i1]*db2[nni1]+da2[i1]*db1[nni1]);
3979 d2zeta22=d2zeta22+coff*(d2a22[i1]*b[nni1]+a[i1]*d2b22[nni1]+
3980 2.0*da2[i1]*db2[nni1]);
3981 d2zeta1r=d2zeta1r+coff*(d2a1r[i1]*b[nni1]+dar[i1]*db1[nni1]+
3982 da1[i1]*dbr[nni1]+a[i1]*d2b1r[nni1]);
3983 d2zeta2r=d2zeta2r+coff*(d2a2r[i1]*b[nni1]+dar[i1]*db2[nni1]+
3984 da2[i1]*dbr[nni1]+a[i1]*d2b2r[nni1]);
3985 deriv2=deriv2+coff*(d2ar2[i1]*b[nni1]+a[i1]*d2br2[nni1]+
3986 2.0*dar[i1]*dbr[nni1]);
3987 }
3988 s=x*0.5;
3989 deriv=0.5*deriv;
3990 deriv2=0.5*deriv2;
3991 dzeta1=0.5*dzeta1;
3992 dzeta2=0.5*dzeta2;
3993 d2zeta11=0.5*d2zeta11;
3994 d2zeta12=0.5*d2zeta12;
3995 d2zeta22=0.5*d2zeta22;
3996 d2zeta1r=0.5*d2zeta1r;
3997 d2zeta2r=0.5*d2zeta2r;
3998
3999 memory->destroy(a);
4000 memory->destroy(b);
4001
4002 return;
4003 }
4004 /* -------------------------------------------------------------------------------
4005 coeffs
4006 ------------------------------------------------------------------------------- */
coeffs(int na,int nb,int k)4007 double PairSMTBQ::coeffs(int na, int nb, int k)
4008 {
4009 // implicit real (a-h,o-z)
4010 // common /fctrl/ fct(30)
4011
4012 int il,je,ia,i,j,ie,l;
4013 double coeffs;
4014
4015 // Statement function
4016 // binm(n,i)=fct(n+1)/(fct(n-i+1)*fct(i+1));
4017
4018 coeffs=0.0;
4019 l=na+nb-k;
4020 ie=min(l,na)+1;
4021 je=min(l,nb);
4022 ia=l-je+1;
4023 for (il = ia; il <= ie; il++) {
4024 i=il-1;
4025 j=l-i; // D'ou vient le i
4026 coeffs=coeffs + binm(na,i)*binm(nb,j)*powint(-1.,j);
4027 }
4028 return coeffs;
4029 }
4030
4031 // ============================================
4032
binm(int n,int i)4033 double PairSMTBQ::binm(int n, int i)
4034 {
4035 return fct[n+1]/(fct[n-i+1]*fct[i+1]);
4036 }
4037
4038 /* ---------------------------------------------------------------------------------
4039 Caintgs
4040 -------------------------------------------------------------------------------- */
caintgs(double x,int k,double * a)4041 void PairSMTBQ::caintgs (double x, int k, double *a)
4042 {
4043 // implicit real (a-h,o-z)
4044 // dimension a(30)
4045 int i;
4046 double cste,rx;
4047
4048 cste=exp(-x);
4049 rx=1.0/x;
4050 a[1]=cste*rx;
4051 for (i = 1; i <= k; i++) {
4052 a[i+1]=(a[i]*static_cast<double>(i)+cste)*rx;
4053 }
4054 return;
4055 }
4056 /* -----------------------------------------------------------------------------------
4057 Cbintgs
4058 ----------------------------------------------------------------------------------- */
cbintgs(double x,int k,double * b)4059 void PairSMTBQ::cbintgs( double x, int k, double *b)
4060 {
4061 // implicit real (a-h,o-z)
4062 /* *******************************************************************
4063 ! Fills array of b-integrals. note that b(i) is b(i-1) in the
4064 ! usual notation
4065 ! for x.gt.3 exponential formula is used
4066 ! for 2.lt.x.le.3 and k.le.10 exponential formula is used
4067 ! for 2.lt.x.le.3 and k.gt.10 15 term series is used
4068 ! for 1.lt.x .e.2 and k.le.7 exponential formula is used
4069 ! for 1.lt.x.le.2 and k.gt.7 12 term series is used
4070 ! for .5.lt.x.le.1 and k.le.5 exponential formula is used
4071 ! for .5.lt.x.le.1 and k.gt.5 7 term series is used
4072 ! for x.le..5 6 term series is used
4073 !******************************************************************* */
4074 // dimension b(30)
4075 // common /fctrl/ fct(30)
4076
4077 int i0,m,last,i;
4078 double absx,expx,expmx,ytrm,y,rx;
4079
4080 i0=0;
4081 absx=fabs(x);
4082
4083 if (absx > 3.0) goto g120;
4084 if (absx > 2.0) goto g20;
4085 if (absx > 1.0) goto g50;
4086 if (absx > 0.5) goto g80;
4087 if (absx > 1.0e-8) goto g110;
4088 goto g170;
4089 g110: last=6;
4090 goto g140;
4091 g80: if (k <= 5) goto g120;
4092 last=7;
4093 goto g140;
4094 g50: if (k <= 7) goto g120;
4095 last=12;
4096 goto g140;
4097 g20: if (k <= 10) goto g120;
4098 last=15;
4099 goto g140;
4100
4101 g120: expx=exp(x);
4102 expmx=1./expx;
4103 rx=1.0/x;
4104 b[1]=(expx-expmx)*rx;
4105 for (i = 1; i <= k ; i++) {
4106 b[i+1]=(static_cast<double>(i)*b[i]+ powint(-1.0,i)*expx-expmx)*rx;
4107 }
4108 goto g190;
4109 //
4110 // Series to calculate b(i)
4111 //
4112 g140: for (i = i0; i <= k ; i++) {
4113 y=0.;
4114 for (m = i0; m <= last; m++) {
4115 ytrm = powint(-x,m-1)*(1. - powint(-1.,m+i+1))/(fct[m+1]*static_cast<double>(m+i+1));
4116 y = y + ytrm*(-x);
4117 }
4118 b[i+1] = y;
4119 }
4120 goto g190;
4121 //
4122 // x extremely small
4123 //
4124 g170: for (i = i0; i <= k; i++) {
4125 b[i+1] = (1.-powint(-1.,i+1))/static_cast<double>(i+1);
4126 }
4127 g190:
4128 return;
4129 }
4130
4131 /* ============================== This is the END... ================================== */
4132