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 #include "neb.h"
16
17 #include "atom.h"
18 #include "domain.h"
19 #include "error.h"
20 #include "finish.h"
21 #include "fix.h"
22 #include "fix_neb.h"
23 #include "math_const.h"
24 #include "memory.h"
25 #include "min.h"
26 #include "modify.h"
27 #include "output.h"
28 #include "thermo.h"
29 #include "timer.h"
30 #include "universe.h"
31 #include "update.h"
32
33 #include <cmath>
34 #include <cstring>
35
36 using namespace LAMMPS_NS;
37 using namespace MathConst;
38
39 #define MAXLINE 256
40 #define CHUNK 1024
41 #define ATTRIBUTE_PERLINE 4
42
43 /* ---------------------------------------------------------------------- */
44
NEB(LAMMPS * lmp)45 NEB::NEB(LAMMPS *lmp) : Command(lmp), all(nullptr), rdist(nullptr) {}
46
47 /* ----------------------------------------------------------------------
48 internal NEB constructor, called from TAD
49 ------------------------------------------------------------------------- */
50
NEB(LAMMPS * lmp,double etol_in,double ftol_in,int n1steps_in,int n2steps_in,int nevery_in,double * buf_init,double * buf_final)51 NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
52 int n2steps_in, int nevery_in, double *buf_init, double *buf_final)
53 : Command(lmp), fp(nullptr), all(nullptr), rdist(nullptr)
54 {
55 double delx,dely,delz;
56
57 etol = etol_in;
58 ftol = ftol_in;
59 n1steps = n1steps_in;
60 n2steps = n2steps_in;
61 nevery = nevery_in;
62 verbose = false;
63
64 // replica info
65
66 nreplica = universe->nworlds;
67 ireplica = universe->iworld;
68 me_universe = universe->me;
69 uworld = universe->uworld;
70 MPI_Comm_rank(world,&me);
71
72 // generate linear interpolate replica
73 double fraction = ireplica/(nreplica-1.0);
74 double **x = atom->x;
75 int nlocal = atom->nlocal;
76
77 int ii = 0;
78 for (int i = 0; i < nlocal; i++) {
79 delx = buf_final[ii] - buf_init[ii];
80 dely = buf_final[ii+1] - buf_init[ii+1];
81 delz = buf_final[ii+2] - buf_init[ii+2];
82 domain->minimum_image(delx,dely,delz);
83 x[i][0] = buf_init[ii] + fraction*delx;
84 x[i][1] = buf_init[ii+1] + fraction*dely;
85 x[i][2] = buf_init[ii+2] + fraction*delz;
86 ii += 3;
87 }
88 }
89
90 /* ---------------------------------------------------------------------- */
91
~NEB()92 NEB::~NEB()
93 {
94 MPI_Comm_free(&roots);
95 memory->destroy(all);
96 delete[] rdist;
97 if (fp) fclose(fp);
98 }
99
100 /* ----------------------------------------------------------------------
101 perform NEB on multiple replicas
102 ------------------------------------------------------------------------- */
103
command(int narg,char ** arg)104 void NEB::command(int narg, char **arg)
105 {
106 if (domain->box_exist == 0)
107 error->all(FLERR,"NEB command before simulation box is defined");
108
109 if (narg < 6) error->universe_all(FLERR,"Illegal NEB command");
110
111 etol = utils::numeric(FLERR,arg[0],false,lmp);
112 ftol = utils::numeric(FLERR,arg[1],false,lmp);
113 n1steps = utils::inumeric(FLERR,arg[2],false,lmp);
114 n2steps = utils::inumeric(FLERR,arg[3],false,lmp);
115 nevery = utils::inumeric(FLERR,arg[4],false,lmp);
116
117 // error checks
118
119 if (etol < 0.0) error->all(FLERR,"Illegal NEB command");
120 if (ftol < 0.0) error->all(FLERR,"Illegal NEB command");
121 if (nevery <= 0) error->universe_all(FLERR,"Illegal NEB command");
122 if (n1steps % nevery || n2steps % nevery)
123 error->universe_all(FLERR,"Illegal NEB command");
124
125 // replica info
126
127 nreplica = universe->nworlds;
128 ireplica = universe->iworld;
129 me_universe = universe->me;
130 uworld = universe->uworld;
131 MPI_Comm_rank(world,&me);
132
133 // error checks
134
135 if (nreplica == 1) error->all(FLERR,"Cannot use NEB with a single replica");
136 if (atom->map_style == Atom::MAP_NONE)
137 error->all(FLERR,"Cannot use NEB unless atom map exists");
138
139 // process file-style setting to setup initial configs for all replicas
140
141 if (strcmp(arg[5],"final") == 0) {
142 if (narg != 7 && narg !=8) error->universe_all(FLERR,"Illegal NEB command");
143 inpfile = arg[6];
144 readfile(inpfile,0);
145 } else if (strcmp(arg[5],"each") == 0) {
146 if (narg != 7 && narg !=8) error->universe_all(FLERR,"Illegal NEB command");
147 inpfile = arg[6];
148 readfile(inpfile,1);
149 } else if (strcmp(arg[5],"none") == 0) {
150 if (narg != 6 && narg !=7) error->universe_all(FLERR,"Illegal NEB command");
151 } else error->universe_all(FLERR,"Illegal NEB command");
152
153 verbose=false;
154 if (strcmp(arg[narg-1],"verbose") == 0) verbose=true;
155 // run the NEB calculation
156
157 run();
158 }
159
160 /* ----------------------------------------------------------------------
161 run NEB on multiple replicas
162 ------------------------------------------------------------------------- */
163
run()164 void NEB::run()
165 {
166 // create MPI communicator for root proc from each world
167
168 int color;
169 if (me == 0) color = 0;
170 else color = 1;
171 MPI_Comm_split(uworld,color,0,&roots);
172
173 int ineb;
174 for (ineb = 0; ineb < modify->nfix; ineb++)
175 if (strcmp(modify->fix[ineb]->style,"neb") == 0) break;
176 if (ineb == modify->nfix) error->all(FLERR,"NEB requires use of fix neb");
177
178 fneb = (FixNEB *) modify->fix[ineb];
179 if (verbose) numall =7;
180 else numall = 4;
181 memory->create(all,nreplica,numall,"neb:all");
182 rdist = new double[nreplica];
183
184 // initialize LAMMPS
185
186 update->whichflag = 2;
187 update->etol = etol;
188 update->ftol = ftol;
189 update->multireplica = 1;
190
191 lmp->init();
192
193 if (update->minimize->searchflag)
194 error->all(FLERR,"NEB requires damped dynamics minimizer");
195
196 // setup regular NEB minimization
197 FILE *uscreen = universe->uscreen;
198 FILE *ulogfile = universe->ulogfile;
199
200 if (me_universe == 0 && uscreen)
201 fprintf(uscreen,"Setting up regular NEB ...\n");
202
203 update->beginstep = update->firststep = update->ntimestep;
204 update->endstep = update->laststep = update->firststep + n1steps;
205 update->nsteps = n1steps;
206 update->max_eval = n1steps;
207 if (update->laststep < 0)
208 error->all(FLERR,"Too many timesteps for NEB");
209
210 update->minimize->setup();
211
212 if (me_universe == 0) {
213 if (uscreen) {
214 if (verbose) {
215 fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
216 "GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
217 "RDN PEN pathangle1 angletangrad1 anglegrad1 gradV1 "
218 "ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
219 "... ReplicaForceN MaxAtomForceN\n");
220 } else {
221 fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
222 "GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
223 "RDN PEN\n");
224 }
225 }
226
227 if (ulogfile) {
228 if (verbose) {
229 fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
230 "GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
231 "RDN PEN pathangle1 angletangrad1 anglegrad1 gradV1 "
232 "ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
233 "... ReplicaForceN MaxAtomForceN\n");
234 } else {
235 fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
236 "GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
237 "RDN PEN\n");
238 }
239 }
240 }
241 print_status();
242
243 // perform regular NEB for n1steps or until replicas converge
244 // retrieve PE values from fix NEB and print every nevery iterations
245 // break out of while loop early if converged
246 // damped dynamic min styles insure all replicas converge together
247
248 timer->init();
249 timer->barrier_start();
250
251 while (update->minimize->niter < n1steps) {
252 update->minimize->run(nevery);
253 print_status();
254 if (update->minimize->stop_condition) break;
255 }
256
257 timer->barrier_stop();
258
259 update->minimize->cleanup();
260
261 Finish finish(lmp);
262 finish.end(1);
263
264 // switch fix NEB to climbing mode
265 // top = replica that becomes hill climber
266
267 double vmax = all[0][0];
268 int top = 0;
269 for (int m = 1; m < nreplica; m++)
270 if (vmax < all[m][0]) {
271 vmax = all[m][0];
272 top = m;
273 }
274
275 // setup climbing NEB minimization
276 // must reinitialize minimizer so it re-creates its fix MINIMIZE
277
278 if (me_universe == 0 && uscreen)
279 fprintf(uscreen,"Setting up climbing ...\n");
280
281 if (me_universe == 0) {
282 if (uscreen)
283 fprintf(uscreen,"Climbing replica = %d\n",top+1);
284 if (ulogfile)
285 fprintf(ulogfile,"Climbing replica = %d\n",top+1);
286 }
287
288 update->beginstep = update->firststep = update->ntimestep;
289 update->endstep = update->laststep = update->firststep + n2steps;
290 update->nsteps = n2steps;
291 update->max_eval = n2steps;
292 if (update->laststep < 0)
293 error->all(FLERR,"Too many timesteps");
294
295 update->minimize->init();
296 fneb->rclimber = top;
297 update->minimize->setup();
298
299 if (me_universe == 0) {
300 if (uscreen) {
301 if (verbose) {
302 fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
303 "GradV0 GradV1 GradVc EBF EBR RDT "
304 "RD1 PE1 RD2 PE2 ... RDN PEN "
305 "pathangle1 angletangrad1 anglegrad1 gradV1 "
306 "ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
307 "... ReplicaForceN MaxAtomForceN\n");
308 } else {
309 fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
310 "GradV0 GradV1 GradVc "
311 "EBF EBR RDT "
312 "RD1 PE1 RD2 PE2 ... RDN PEN\n");
313 }
314 }
315 if (ulogfile) {
316 if (verbose) {
317 fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
318 "GradV0 GradV1 GradVc EBF EBR RDT "
319 "RD1 PE1 RD2 PE2 ... RDN PEN "
320 "pathangle1 angletangrad1 anglegrad1 gradV1 "
321 "ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
322 "... ReplicaForceN MaxAtomForceN\n");
323 } else {
324 fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
325 "GradV0 GradV1 GradVc "
326 "EBF EBR RDT "
327 "RD1 PE1 RD2 PE2 ... RDN PEN\n");
328 }
329 }
330 }
331 print_status();
332
333 // perform climbing NEB for n2steps or until replicas converge
334 // retrieve PE values from fix NEB and print every nevery iterations
335 // break induced if converged
336 // damped dynamic min styles insure all replicas converge together
337
338 timer->init();
339 timer->barrier_start();
340
341 while (update->minimize->niter < n2steps) {
342 update->minimize->run(nevery);
343 print_status();
344 if (update->minimize->stop_condition) break;
345 }
346
347 timer->barrier_stop();
348
349 update->minimize->cleanup();
350
351 finish.end(1);
352
353 update->whichflag = 0;
354 update->multireplica = 0;
355 update->firststep = update->laststep = 0;
356 update->beginstep = update->endstep = 0;
357 }
358
359 /* ----------------------------------------------------------------------
360 read initial config atom coords from file
361 flag = 0
362 only first replica opens file and reads it
363 first replica bcasts lines to all replicas
364 final replica stores coords
365 intermediate replicas interpolate from coords
366 new coord = replica fraction between current and final state
367 initial replica does nothing
368 flag = 1
369 each replica (except first) opens file and reads it
370 each replica stores coords
371 initial replica does nothing
372 ------------------------------------------------------------------------- */
373
readfile(char * file,int flag)374 void NEB::readfile(char *file, int flag)
375 {
376 int i,j,m,nchunk,eofflag,nlines;
377 tagint tag;
378 char *eof,*start,*next,*buf;
379 char line[MAXLINE];
380 double xx,yy,zz,delx,dely,delz;
381
382 if (me_universe == 0 && universe->uscreen)
383 fprintf(universe->uscreen,"Reading NEB coordinate file(s) ...\n");
384
385 // flag = 0, universe root reads header of file, bcast to universe
386 // flag = 1, each replica's root reads header of file, bcast to world
387 // but explicitly skip first replica
388
389 if (flag == 0) {
390 if (me_universe == 0) {
391 open(file);
392 while (1) {
393 eof = fgets(line,MAXLINE,fp);
394 if (eof == nullptr) error->one(FLERR,"Unexpected end of NEB file");
395 start = &line[strspn(line," \t\n\v\f\r")];
396 if (*start != '\0' && *start != '#') break;
397 }
398 int rv = sscanf(line,"%d",&nlines);
399 if (rv != 1) nlines = -1;
400 }
401 MPI_Bcast(&nlines,1,MPI_INT,0,uworld);
402 if (nlines < 0)
403 error->universe_all(FLERR,"Incorrectly formatted NEB file");
404 } else {
405 if (me == 0) {
406 if (ireplica) {
407 open(file);
408 while (1) {
409 eof = fgets(line,MAXLINE,fp);
410 if (eof == nullptr) error->one(FLERR,"Unexpected end of NEB file");
411 start = &line[strspn(line," \t\n\v\f\r")];
412 if (*start != '\0' && *start != '#') break;
413 }
414 int rv = sscanf(line,"%d",&nlines);
415 if (rv != 1) nlines = -1;
416 } else nlines = 0;
417 }
418 MPI_Bcast(&nlines,1,MPI_INT,0,world);
419 if (nlines < 0)
420 error->all(FLERR,"Incorrectly formatted NEB file");
421 }
422
423 char *buffer = new char[CHUNK*MAXLINE];
424 char **values = new char*[ATTRIBUTE_PERLINE];
425
426 double fraction = ireplica/(nreplica-1.0);
427
428 double **x = atom->x;
429 int nlocal = atom->nlocal;
430
431 // loop over chunks of lines read from file
432 // two versions of read_lines_from_file() for world vs universe bcast
433 // count # of atom coords changed so can check for invalid atom IDs in file
434
435 int ncount = 0;
436
437 int nread = 0;
438 while (nread < nlines) {
439 nchunk = MIN(nlines-nread,CHUNK);
440 if (flag == 0)
441 eofflag = utils::read_lines_from_file(fp,nchunk,MAXLINE,buffer,
442 universe->me,universe->uworld);
443 else
444 eofflag = utils::read_lines_from_file(fp,nchunk,MAXLINE,buffer,me,world);
445 if (eofflag) error->all(FLERR,"Unexpected end of NEB file");
446
447 buf = buffer;
448 next = strchr(buf,'\n');
449 *next = '\0';
450 int nwords = utils::count_words(utils::trim_comment(buf));
451 *next = '\n';
452
453 if (nwords != ATTRIBUTE_PERLINE)
454 error->all(FLERR,"Incorrect atom format in NEB file");
455
456 // loop over lines of atom coords
457 // tokenize the line into values
458
459 for (i = 0; i < nchunk; i++) {
460 next = strchr(buf,'\n');
461
462 values[0] = strtok(buf," \t\n\r\f");
463 for (j = 1; j < nwords; j++)
464 values[j] = strtok(nullptr," \t\n\r\f");
465
466 // adjust atom coord based on replica fraction
467 // for flag = 0, interpolate for intermediate and final replicas
468 // for flag = 1, replace existing coord with new coord
469 // ignore image flags of replica x
470 // displacement from first replica is via minimum image convention
471 // if x of some replica is across periodic boundary:
472 // new x may be outside box
473 // will be remapped back into box when simulation starts
474 // its image flags will then be adjusted
475
476 tag = ATOTAGINT(values[0]);
477 m = atom->map(tag);
478 if (m >= 0 && m < nlocal) {
479 ncount++;
480 xx = atof(values[1]);
481 yy = atof(values[2]);
482 zz = atof(values[3]);
483
484 delx = xx - x[m][0];
485 dely = yy - x[m][1];
486 delz = zz - x[m][2];
487
488 domain->minimum_image(delx,dely,delz);
489
490 if (flag == 0) {
491 x[m][0] += fraction*delx;
492 x[m][1] += fraction*dely;
493 x[m][2] += fraction*delz;
494 } else {
495 x[m][0] += delx;
496 x[m][1] += dely;
497 x[m][2] += delz;
498 }
499 }
500
501 buf = next + 1;
502 }
503
504 nread += nchunk;
505 }
506
507 // check that all atom IDs in file were found by a proc
508
509 if (flag == 0) {
510 int ntotal;
511 MPI_Allreduce(&ncount,&ntotal,1,MPI_INT,MPI_SUM,uworld);
512 if (ntotal != nreplica*nlines)
513 error->universe_all(FLERR,"Invalid atom IDs in NEB file");
514 } else {
515 int ntotal;
516 MPI_Allreduce(&ncount,&ntotal,1,MPI_INT,MPI_SUM,world);
517 if (ntotal != nlines)
518 error->all(FLERR,"Invalid atom IDs in NEB file");
519 }
520
521 // clean up
522
523 delete [] buffer;
524 delete [] values;
525
526 if (flag == 0) {
527 if (me_universe == 0) {
528 if (compressed) pclose(fp);
529 else fclose(fp);
530 }
531 } else {
532 if (me == 0 && ireplica) {
533 if (compressed) pclose(fp);
534 else fclose(fp);
535 }
536 }
537 fp = nullptr;
538 }
539
540 /* ----------------------------------------------------------------------
541 universe proc 0 opens NEB data file
542 test if gzipped
543 ------------------------------------------------------------------------- */
544
open(char * file)545 void NEB::open(char *file)
546 {
547 compressed = 0;
548 char *suffix = file + strlen(file) - 3;
549 if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1;
550 if (!compressed) fp = fopen(file,"r");
551 else {
552 #ifdef LAMMPS_GZIP
553 auto gunzip = std::string("gzip -c -d ") + file;
554 #ifdef _WIN32
555 fp = _popen(gunzip.c_str(),"rb");
556 #else
557 fp = popen(gunzip.c_str(),"r");
558 #endif
559
560 #else
561 error->one(FLERR,"Cannot open gzipped file");
562 #endif
563 }
564
565 if (fp == nullptr)
566 error->one(FLERR,"Cannot open file {}: {}",file,utils::getsyserror());
567 }
568
569 /* ----------------------------------------------------------------------
570 query fix NEB for info on each replica
571 universe proc 0 prints current NEB status
572 ------------------------------------------------------------------------- */
573
print_status()574 void NEB::print_status()
575 {
576 double fnorm2 = sqrt(update->minimize->fnorm_sqr());
577 double fmaxreplica;
578 MPI_Allreduce(&fnorm2,&fmaxreplica,1,MPI_DOUBLE,MPI_MAX,roots);
579 double fnorminf = update->minimize->fnorm_inf();
580 double fmaxatom;
581 MPI_Allreduce(&fnorminf,&fmaxatom,1,MPI_DOUBLE,MPI_MAX,roots);
582
583 if (verbose) {
584 freplica = new double[nreplica];
585 MPI_Allgather(&fnorm2,1,MPI_DOUBLE,&freplica[0],1,MPI_DOUBLE,roots);
586 fmaxatomInRepl = new double[nreplica];
587 MPI_Allgather(&fnorminf,1,MPI_DOUBLE,&fmaxatomInRepl[0],1,MPI_DOUBLE,roots);
588 }
589
590 double one[7];
591 one[0] = fneb->veng;
592 one[1] = fneb->plen;
593 one[2] = fneb->nlen;
594 one[3] = fneb->gradlen;
595
596 if (verbose) {
597 one[4] = fneb->dotpath;
598 one[5] = fneb->dottangrad;
599 one[6] = fneb->dotgrad;
600 }
601
602 if (output->thermo->normflag) one[0] /= atom->natoms;
603 if (me == 0)
604 MPI_Allgather(one,numall,MPI_DOUBLE,&all[0][0],numall,MPI_DOUBLE,roots);
605 MPI_Bcast(&all[0][0],numall*nreplica,MPI_DOUBLE,0,world);
606
607 rdist[0] = 0.0;
608 for (int i = 1; i < nreplica; i++)
609 rdist[i] = rdist[i-1] + all[i][1];
610 double endpt = rdist[nreplica-1] = rdist[nreplica-2] + all[nreplica-2][2];
611 for (int i = 1; i < nreplica; i++)
612 rdist[i] /= endpt;
613
614 // look up GradV for the initial, final, and climbing replicas
615 // these are identical to fnorm2, but to be safe we
616 // take them straight from fix_neb
617
618 double gradvnorm0, gradvnorm1, gradvnormc;
619
620 int irep;
621 irep = 0;
622 gradvnorm0 = all[irep][3];
623 irep = nreplica-1;
624 gradvnorm1 = all[irep][3];
625 irep = fneb->rclimber;
626 if (irep > -1) {
627 gradvnormc = all[irep][3];
628 ebf = all[irep][0]-all[0][0];
629 ebr = all[irep][0]-all[nreplica-1][0];
630 } else {
631 double vmax = all[0][0];
632 int top = 0;
633 for (int m = 1; m < nreplica; m++)
634 if (vmax < all[m][0]) {
635 vmax = all[m][0];
636 top = m;
637 }
638 irep = top;
639 gradvnormc = all[irep][3];
640 ebf = all[irep][0]-all[0][0];
641 ebr = all[irep][0]-all[nreplica-1][0];
642 }
643
644 if (me_universe == 0) {
645 const double todeg=180.0/MY_PI;
646 FILE *uscreen = universe->uscreen;
647 FILE *ulogfile = universe->ulogfile;
648 if (uscreen) {
649 fprintf(uscreen,BIGINT_FORMAT " %12.8g %12.8g ",
650 update->ntimestep,fmaxreplica,fmaxatom);
651 fprintf(uscreen,"%12.8g %12.8g %12.8g ",
652 gradvnorm0,gradvnorm1,gradvnormc);
653 fprintf(uscreen,"%12.8g %12.8g %12.8g ",ebf,ebr,endpt);
654 for (int i = 0; i < nreplica; i++)
655 fprintf(uscreen,"%12.8g %12.8g ",rdist[i],all[i][0]);
656 if (verbose) {
657 fprintf(uscreen,"%12.5g %12.5g %12.5g %12.5g %12.5g %12.5g",
658 NAN,180-acos(all[0][5])*todeg,180-acos(all[0][6])*todeg,
659 all[0][3],freplica[0],fmaxatomInRepl[0]);
660 for (int i = 1; i < nreplica-1; i++)
661 fprintf(uscreen,"%12.5g %12.5g %12.5g %12.5g %12.5g %12.5g",
662 180-acos(all[i][4])*todeg,180-acos(all[i][5])*todeg,
663 180-acos(all[i][6])*todeg,all[i][3],freplica[i],
664 fmaxatomInRepl[i]);
665 fprintf(uscreen,"%12.5g %12.5g %12.5g %12.5g %12.5g %12.5g",
666 NAN,180-acos(all[nreplica-1][5])*todeg,NAN,all[nreplica-1][3],
667 freplica[nreplica-1],fmaxatomInRepl[nreplica-1]);
668 }
669 fprintf(uscreen,"\n");
670 }
671
672 if (ulogfile) {
673 fprintf(ulogfile,BIGINT_FORMAT " %12.8g %12.8g ",
674 update->ntimestep,fmaxreplica,fmaxatom);
675 fprintf(ulogfile,"%12.8g %12.8g %12.8g ",
676 gradvnorm0,gradvnorm1,gradvnormc);
677 fprintf(ulogfile,"%12.8g %12.8g %12.8g ",ebf,ebr,endpt);
678 for (int i = 0; i < nreplica; i++)
679 fprintf(ulogfile,"%12.8g %12.8g ",rdist[i],all[i][0]);
680 if (verbose) {
681 fprintf(ulogfile,"%12.5g %12.5g %12.5g %12.5g %12.5g %12.5g",
682 NAN,180-acos(all[0][5])*todeg,180-acos(all[0][6])*todeg,
683 all[0][3],freplica[0],fmaxatomInRepl[0]);
684 for (int i = 1; i < nreplica-1; i++)
685 fprintf(ulogfile,"%12.5g %12.5g %12.5g %12.5g %12.5g %12.5g",
686 180-acos(all[i][4])*todeg,180-acos(all[i][5])*todeg,
687 180-acos(all[i][6])*todeg,all[i][3],freplica[i],
688 fmaxatomInRepl[i]);
689 fprintf(ulogfile,"%12.5g %12.5g %12.5g %12.5g %12.5g %12.5g",
690 NAN,180-acos(all[nreplica-1][5])*todeg,NAN,all[nreplica-1][3],
691 freplica[nreplica-1],fmaxatomInRepl[nreplica-1]);
692 }
693 fprintf(ulogfile,"\n");
694 fflush(ulogfile);
695 }
696 }
697 }
698