← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to src/USER-COLVARS/fix_colvars.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2013-11-20 22:41:36 UTC
  • mfrom: (1.2.2)
  • Revision ID: package-import@ubuntu.com-20131120224136-tzx7leh606fqnckm
Tags: 0~20131119.git7162cf0-1
* [e65b919] Imported Upstream version 0~20131119.git7162cf0
* [f7bddd4] Fix some problems, introduced by upstream recently.
* [3616dfc] Use wrap-and-sort script.
* [7e92030] Ignore quilt dir

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/USER-COLVARS/fix_colvars.cpp
10
10
 
11
11
   See the README file in the top-level LAMMPS directory.
12
12
------------------------------------------------------------------------- */
 
13
 
13
14
/* ----------------------------------------------------------------------
14
 
   Contributing author:  Axel Kohlmeyer (TempleU)
 
15
   Contributing author:  Axel Kohlmeyer (Temple U)
15
16
------------------------------------------------------------------------- */
16
17
 
17
18
#include <math.h>
18
19
#include <stdio.h>
19
20
#include <stdlib.h>
20
21
#include <string.h>
21
 
 
22
22
#include <errno.h>
23
23
 
24
24
#include "fix_colvars.h"
31
31
#include "modify.h"
32
32
#include "respa.h"
33
33
#include "update.h"
 
34
#include "citeme.h"
34
35
 
35
36
#include "colvarproxy_lammps.h"
36
37
 
 
38
static const char colvars_pub[] =
 
39
  "fix colvars command:\n\n"
 
40
  "@Article{fiorin13,\n"
 
41
  " author =  {G.~Fiorin and M.{\\,}L.~Klein and J.~H{\\'e}nin},\n"
 
42
  " title =   {Using collective variables to drive molecular"
 
43
  " dynamics simulations},\n"
 
44
  " journal = {Mol.~Phys.},\n"
 
45
  " year =    2013,\n"
 
46
  " note =    {doi: 10.1080/00268976.2013.813594}\n"
 
47
  "}\n\n";
 
48
 
37
49
/* re-usable integer hash table code with static linkage. */
38
50
 
39
51
/** hash table top level data structure */
286
298
  seed    <integer>         (seed for RNG, defaults to '1966')
287
299
  tstat   <fix label>       (label of thermostatting fix)
288
300
 
289
 
 TODO: add (optional) arguments for RNG seed, temperature compute
290
301
 ***************************************************************/
291
302
FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) :
292
303
  Fix(lmp, narg, arg)
294
305
  if (narg < 4)
295
306
    error->all(FLERR,"Illegal fix colvars command: too few arguments");
296
307
 
297
 
  if (atom->tag_enable == 0)
298
 
    error->all(FLERR,"Cannot use fix colvars without atom IDs defined");
299
 
 
300
 
  if (atom->rmass_flag)
301
 
    error->all(FLERR,"Cannot use fix colvars for atoms with rmass attribute");
302
 
 
303
308
  if (instances > 0)
304
 
    error->all(FLERR,"Only one fix colvars can be active at a time");
 
309
    error->all(FLERR,"Only one colvars fix can be active at a time");
305
310
  ++instances;
306
311
 
307
312
  scalar_flag = 1;
308
313
  global_freq = 1;
309
314
  nevery = 1;
310
315
  extscalar = 1;
 
316
  restart_global = 1;
311
317
 
312
318
  me = comm->me;
313
319
 
314
320
  conf_file = strdup(arg[3]);
315
321
  rng_seed = 1966;
316
 
  unwrap_flag = 0;
 
322
  unwrap_flag = 1;
317
323
 
318
324
  inp_name = NULL;
319
325
  out_name = NULL;
321
327
 
322
328
  /* parse optional arguments */
323
329
  int argsdone = 4;
324
 
  while (argsdone+1 < narg) {
 
330
  while (argsdone < narg) {
 
331
    // we have keyword/value pairs. check if value is missing
 
332
    if (argsdone+1 == narg)
 
333
      error->all(FLERR,"Missing argument to keyword");
 
334
 
325
335
    if (0 == strcmp(arg[argsdone], "input")) {
326
336
      inp_name = strdup(arg[argsdone+1]);
327
337
    } else if (0 == strcmp(arg[argsdone], "output")) {
328
338
      out_name = strdup(arg[argsdone+1]);
329
339
    } else if (0 == strcmp(arg[argsdone], "seed")) {
330
 
      rng_seed = atoi(arg[argsdone+1]);
 
340
      rng_seed = force->inumeric(FLERR,arg[argsdone+1]);
331
341
    } else if (0 == strcmp(arg[argsdone], "unwrap")) {
332
342
      if (0 == strcmp(arg[argsdone+1], "yes")) {
333
343
        unwrap_flag = 1;
350
360
  tstat_id = -1;
351
361
  energy = 0.0;
352
362
  nlevels_respa = 0;
 
363
  init_flag = 0;
353
364
  num_coords = 0;
354
365
  coords = forces = oforce = NULL;
355
366
  comm_buf = NULL;
359
370
 
360
371
  /* storage required to communicate a single coordinate or force. */
361
372
  size_one = sizeof(struct commdata);
 
373
 
 
374
  if (lmp->citeme) lmp->citeme->add(colvars_pub);
362
375
}
363
376
 
364
377
/*********************************
370
383
  memory->sfree(inp_name);
371
384
  memory->sfree(out_name);
372
385
  memory->sfree(tmp_name);
373
 
  deallocate();
374
 
  --instances;
375
 
}
376
 
 
377
 
/* ---------------------------------------------------------------------- */
378
 
 
379
 
void FixColvars::deallocate()
380
 
{
381
 
  memory->destroy(comm_buf);
 
386
  memory->sfree(comm_buf);
382
387
 
383
388
  if (proxy) {
384
389
    delete proxy;
387
392
    delete hashtable;
388
393
  }
389
394
 
390
 
  proxy = NULL;
391
 
  idmap = NULL;
392
 
  coords = NULL;
393
 
  forces = NULL;
394
 
  oforce = NULL;
395
 
  comm_buf = NULL;
396
 
}
397
 
 
398
 
/* ---------------------------------------------------------------------- */
399
 
 
400
 
void FixColvars::post_run()
401
 
{
402
 
  deallocate();
403
 
  memory->sfree(inp_name);
404
 
  inp_name = strdup(out_name);
 
395
  --instances;
405
396
}
406
397
 
407
398
/* ---------------------------------------------------------------------- */
413
404
  mask |= MIN_POST_FORCE;
414
405
  mask |= POST_FORCE;
415
406
  mask |= POST_FORCE_RESPA;
416
 
  mask |= POST_RUN;
417
407
  mask |= END_OF_STEP;
418
408
  return mask;
419
409
}
420
410
 
421
411
/* ---------------------------------------------------------------------- */
422
412
 
423
 
// initial setup of colvars run.
 
413
// initial checks for colvars run.
424
414
 
425
415
void FixColvars::init()
426
416
{
 
417
  if (atom->tag_enable == 0)
 
418
    error->all(FLERR,"Cannot use fix colvars without atom IDs");
 
419
 
 
420
  if (atom->map_style == 0)
 
421
    error->all(FLERR,"Fix colvars requires an atom map");
 
422
 
 
423
  if ((me == 0) && (update->whichflag == 2))
 
424
    error->warning(FLERR,"Using fix colvars with minimization");
 
425
 
427
426
  if (strstr(update->integrate_style,"respa"))
428
427
    nlevels_respa = ((Respa *) update->integrate)->nlevels;
429
 
 
430
 
 
431
 
  int i,nme,tmp,ndata;
432
 
 
433
 
  MPI_Status status;
434
 
  MPI_Request request;
435
 
 
436
 
  // collect a list of atom type by atom id for the entire system.
437
 
  // the colvar module requires this information to set masses. :-(
438
 
 
439
 
  int *typemap,*type_buf;
440
 
  int nlocal_max,tag_max,max;
441
 
  const int * const tag  = atom->tag;
442
 
  const int * const type = atom->type;
443
 
  int nlocal = atom->nlocal;
444
 
 
445
 
  max=0;
446
 
  for (i = 0; i < nlocal; i++) max = MAX(max,tag[i]);
447
 
  MPI_Allreduce(&max,&tag_max,1,MPI_INT,MPI_MAX,world);
448
 
  MPI_Allreduce(&nlocal,&nlocal_max,1,MPI_INT,MPI_MAX,world);
449
 
 
450
 
  if (me == 0) {
451
 
    typemap = new int[tag_max+1];
452
 
    memset(typemap,0,sizeof(int)*tag_max);
453
 
  }
454
 
  type_buf = new int[2*nlocal_max];
455
 
 
456
 
  if (me == 0) {
457
 
    for (i=0; i<nlocal; ++i)
458
 
      typemap[tag[i]] = type[i];
459
 
 
460
 
    // loop over procs to receive and apply remote data
461
 
 
462
 
    for (i=1; i < comm->nprocs; ++i) {
463
 
      MPI_Irecv(type_buf, 2*nlocal_max, MPI_INT, i, 0, world, &request);
464
 
      MPI_Send(&tmp, 0, MPI_INT, i, 0, world);
465
 
      MPI_Wait(&request, &status);
466
 
      MPI_Get_count(&status, MPI_INT, &ndata);
467
 
 
468
 
      for (int k=0; k<ndata; k+=2)
469
 
        typemap[type_buf[k]] = type_buf[k+1];
470
 
    }
471
 
  } else { // me != 0
472
 
 
473
 
    // copy tag/type data into communication buffer
474
 
 
475
 
    nme = 0;
476
 
    for (i=0; i<nlocal; ++i) {
477
 
      type_buf[nme] = tag[i];
478
 
      type_buf[nme+1] = type[i];
479
 
      nme +=2;
480
 
    }
481
 
    /* blocking receive to wait until it is our turn to send data. */
482
 
    MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status);
483
 
    MPI_Rsend(type_buf, nme, MPI_INT, 0, 0, world);
484
 
  }
485
 
 
486
 
  // now create and initialize the colvar proxy
487
 
 
488
 
  if (me == 0) {
 
428
}
 
429
 
 
430
 
 
431
/* ---------------------------------------------------------------------- */
 
432
 
 
433
void FixColvars::one_time_init()
 
434
{
 
435
  int i,tmp;
 
436
 
 
437
  if (init_flag) return;
 
438
  init_flag = 1;
 
439
 
 
440
   // create and initialize the colvars proxy
 
441
 
 
442
  if (me == 0) {
 
443
    if (screen) fputs("colvars: Creating proxy instance\n",screen);
 
444
    if (logfile) fputs("colvars: Creating proxy instance\n",logfile);
489
445
 
490
446
    if (inp_name) {
491
447
      if (strcmp(inp_name,"NULL") == 0) {
494
450
      }
495
451
    }
496
452
 
 
453
    // try to determine thermostat target temperature
497
454
    double t_target = 0.0;
498
455
    if (tmp_name) {
499
 
      if (strcmp(tmp_name,"NULL") == 0)
 
456
      if (strcmp(tmp_name,"NULL") == 0) {
500
457
        tstat_id = -1;
501
 
      else {
 
458
      } else {
502
459
        tstat_id = modify->find_fix(tmp_name);
503
460
        if (tstat_id < 0) error->one(FLERR,"Could not find tstat fix ID");
504
461
        double *tt = (double*)modify->fix[tstat_id]->extract("t_target",tmp);
506
463
      }
507
464
    }
508
465
 
509
 
    proxy = new colvarproxy_lammps(lmp,conf_file,inp_name,out_name,
510
 
                                   rng_seed,t_target,typemap);
 
466
    proxy = new colvarproxy_lammps(lmp,inp_name,out_name,rng_seed,t_target);
 
467
    proxy->init(conf_file);
511
468
    coords = proxy->get_coords();
512
469
    forces = proxy->get_forces();
513
470
    oforce = proxy->get_oforce();
534
491
    }
535
492
  }
536
493
  MPI_Bcast(taglist, num_coords, MPI_INT, 0, world);
 
494
}
 
495
 
 
496
/* ---------------------------------------------------------------------- */
 
497
 
 
498
void FixColvars::setup(int vflag)
 
499
{
 
500
  const int * const tag  = atom->tag;
 
501
  const int * const type = atom->type;
 
502
  int i,nme,tmp,ndata;
 
503
  int nlocal = atom->nlocal;
 
504
 
 
505
  MPI_Status status;
 
506
  MPI_Request request;
 
507
 
 
508
  one_time_init();
537
509
 
538
510
  // determine size of comm buffer
539
511
  nme=0;
584
556
          cd[i].y = x[k][1];
585
557
          cd[i].z = x[k][2];
586
558
        }
 
559
        if (atom->rmass_flag) {
 
560
          cd[i].m = atom->rmass[k];
 
561
        } else {
 
562
          cd[i].m = atom->mass[type[k]];
 
563
        }
587
564
      }
588
565
    }
589
566
 
605
582
          cd[j].x = comm_buf[k].x;
606
583
          cd[j].y = comm_buf[k].y;
607
584
          cd[j].z = comm_buf[k].z;
 
585
          cd[j].m = comm_buf[k].m;
608
586
          of[j].x = of[j].y = of[j].z = 0.0;
609
587
        }
610
588
      }
635
613
          comm_buf[nme].z = x[k][2];
636
614
        }
637
615
 
 
616
        if (atom->rmass_flag) {
 
617
          comm_buf[nme].m = atom->rmass[k];
 
618
        } else {
 
619
          comm_buf[nme].m = atom->mass[type[k]];
 
620
        }
 
621
 
638
622
        ++nme;
639
623
      }
640
624
    }
643
627
    MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world);
644
628
  }
645
629
 
646
 
  // clear temporary storage
647
 
  if (me == 0) delete typemap;
648
 
  delete type_buf;
649
 
}
650
 
 
651
 
/* ---------------------------------------------------------------------- */
652
 
 
653
 
void FixColvars::setup(int vflag)
654
 
{
655
 
  if (strstr(update->integrate_style,"verlet"))
 
630
  // run pre-run setup in colvarproxy
 
631
  if (me == 0)
 
632
    proxy->setup();
 
633
 
 
634
  // initialize forces
 
635
  if (strstr(update->integrate_style,"verlet") || (update->whichflag == 2))
656
636
    post_force(vflag);
657
 
  else
658
 
    post_force_respa(vflag,0,0);
 
637
  else {
 
638
    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
 
639
    post_force_respa(vflag,nlevels_respa-1,0);
 
640
    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
 
641
  }
659
642
}
660
643
 
661
644
/* ---------------------------------------------------------------------- */
666
649
  // some housekeeping: update status of the proxy as needed.
667
650
  if (me == 0) {
668
651
    if (proxy->want_exit())
669
 
      error->one(FLERR,"Run halted on request from colvars module.\n");
 
652
      error->one(FLERR,"Run aborted on request from colvars module.\n");
670
653
 
671
654
    if (tstat_id < 0) {
672
655
      proxy->set_temperature(0.0);
922
905
 
923
906
/* ---------------------------------------------------------------------- */
924
907
 
 
908
void FixColvars::write_restart(FILE *fp)
 
909
{
 
910
  if (me == 0) {
 
911
    std::string rest_text("");
 
912
    proxy->serialize_status(rest_text);
 
913
    const char *cvm_state = rest_text.c_str();
 
914
    int len = strlen(cvm_state) + 1; // need to include terminating NULL byte.
 
915
    fwrite(&len,sizeof(int),1,fp);
 
916
    fwrite(cvm_state,1,len,fp);
 
917
  }
 
918
}
 
919
 
 
920
/* ---------------------------------------------------------------------- */
 
921
 
 
922
void FixColvars::restart(char *buf)
 
923
{
 
924
  one_time_init();
 
925
 
 
926
  if (me == 0) {
 
927
    std::string rest_text(buf);
 
928
    proxy->deserialize_status(rest_text);
 
929
  }
 
930
}
 
931
 
 
932
/* ---------------------------------------------------------------------- */
 
933
 
925
934
double FixColvars::compute_scalar()
926
935
{
927
936
  return energy;