← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/USER-FEP/fix_adapt_fep.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/USER-FEP/fix_adapt_fep.cpp
20
20
#include "stdlib.h"
21
21
#include "fix_adapt_fep.h"
22
22
#include "atom.h"
23
 
#include "comm.h"
24
23
#include "update.h"
 
24
#include "group.h"
25
25
#include "modify.h"
26
26
#include "force.h"
27
27
#include "pair.h"
28
28
#include "pair_hybrid.h"
29
29
#include "kspace.h"
 
30
#include "fix_store.h"
30
31
#include "input.h"
31
32
#include "variable.h"
 
33
#include "respa.h"
32
34
#include "math_const.h"
33
35
#include "memory.h"
34
36
#include "error.h"
40
42
enum{PAIR,KSPACE,ATOM};
41
43
enum{DIAMETER,CHARGE};
42
44
 
43
 
#undef ADAPT_DEBUG
44
 
 
45
45
/* ---------------------------------------------------------------------- */
46
46
 
47
47
FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
51
51
  nevery = force->inumeric(FLERR,arg[3]);
52
52
  if (nevery < 0) error->all(FLERR,"Illegal fix adapt/fep command");
53
53
 
 
54
  dynamic_group_allow = 1;
 
55
  create_attribute = 1;
 
56
 
54
57
  // count # of adaptations
55
58
 
56
59
  nadapt = 0;
169
172
  for (int m = 0; m < nadapt; m++)
170
173
    if (adapt[m].which == PAIR)
171
174
      memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
 
175
 
 
176
  id_fix_diam = id_fix_chg = NULL;
172
177
}
173
178
 
174
179
/* ---------------------------------------------------------------------- */
184
189
    }
185
190
  }
186
191
  delete [] adapt;
 
192
 
 
193
  // check nfix in case all fixes have already been deleted
 
194
 
 
195
  if (id_fix_diam && modify->nfix) modify->delete_fix(id_fix_diam);
 
196
  if (id_fix_chg && modify->nfix) modify->delete_fix(id_fix_chg);
 
197
  delete [] id_fix_diam;
 
198
  delete [] id_fix_chg;
187
199
}
188
200
 
189
201
/* ---------------------------------------------------------------------- */
193
205
  int mask = 0;
194
206
  mask |= PRE_FORCE;
195
207
  mask |= POST_RUN;
 
208
  mask |= PRE_FORCE_RESPA;
196
209
  return mask;
197
210
}
198
211
 
 
212
/* ----------------------------------------------------------------------
 
213
   if need to restore per-atom quantities, create new fix STORE styles
 
214
------------------------------------------------------------------------- */
 
215
 
 
216
void FixAdaptFEP::post_constructor()
 
217
{
 
218
  if (!resetflag) return;
 
219
  if (!diamflag && !chgflag) return;
 
220
 
 
221
  // new id = fix-ID + FIX_STORE_ATTRIBUTE
 
222
  // new fix group = group for this fix
 
223
 
 
224
  id_fix_diam = NULL;
 
225
  id_fix_chg = NULL;
 
226
 
 
227
  char **newarg = new char*[5];
 
228
  newarg[1] = group->names[igroup];
 
229
  newarg[2] = (char *) "STORE";
 
230
  newarg[3] = (char *) "1";
 
231
  newarg[4] = (char *) "1";
 
232
 
 
233
  if (diamflag) {
 
234
    int n = strlen(id) + strlen("_FIX_STORE_DIAM") + 1;
 
235
    id_fix_diam = new char[n];
 
236
    strcpy(id_fix_diam,id);
 
237
    strcat(id_fix_diam,"_FIX_STORE_DIAM");
 
238
    newarg[0] = id_fix_diam;
 
239
    modify->add_fix(5,newarg);
 
240
    fix_diam = (FixStore *) modify->fix[modify->nfix-1];
 
241
 
 
242
    if (fix_diam->restart_reset) fix_diam->restart_reset = 0;
 
243
    else {
 
244
      double *vec = fix_diam->vstore;
 
245
      double *radius = atom->radius;
 
246
      int *mask = atom->mask;
 
247
      int nlocal = atom->nlocal;
 
248
 
 
249
      for (int i = 0; i < nlocal; i++) {
 
250
        if (mask[i] & groupbit) vec[i] = radius[i];
 
251
        else vec[i] = 0.0;
 
252
      }
 
253
    }
 
254
  }
 
255
 
 
256
  if (chgflag) {
 
257
    int n = strlen(id) + strlen("_FIX_STORE_CHG") + 1;
 
258
    id_fix_chg = new char[n];
 
259
    strcpy(id_fix_chg,id);
 
260
    strcat(id_fix_chg,"_FIX_STORE_CHG");
 
261
    newarg[0] = id_fix_chg;
 
262
    modify->add_fix(5,newarg);
 
263
    fix_chg = (FixStore *) modify->fix[modify->nfix-1];
 
264
 
 
265
    if (fix_chg->restart_reset) fix_chg->restart_reset = 0;
 
266
    else {
 
267
      double *vec = fix_chg->vstore;
 
268
      double *q = atom->q;
 
269
      int *mask = atom->mask;
 
270
      int nlocal = atom->nlocal;
 
271
 
 
272
      for (int i = 0; i < nlocal; i++) {
 
273
        if (mask[i] & groupbit) vec[i] = q[i];
 
274
        else vec[i] = 0.0;
 
275
      }
 
276
    }
 
277
  }
 
278
 
 
279
  delete [] newarg;
 
280
}
 
281
 
199
282
/* ---------------------------------------------------------------------- */
200
283
 
201
284
void FixAdaptFEP::init()
202
285
{
203
286
  int i,j;
204
287
 
 
288
  // allow a dynamic group only if ATOM attribute not used
 
289
 
 
290
  if (group->dynamic[igroup])
 
291
    for (int i = 0; i < nadapt; i++)
 
292
      if (adapt[i].which == ATOM) 
 
293
        error->all(FLERR,"Cannot use dynamic group with fix adapt/fep atom");
 
294
 
205
295
  // setup and error checks
206
296
 
207
297
  anypair = 0;
217
307
 
218
308
    if (ad->which == PAIR) {
219
309
      anypair = 1;
 
310
      Pair *pair = NULL;
220
311
 
221
 
      Pair *pair = force->pair_match(ad->pstyle,1);
222
 
      if (pair == NULL) error->all(FLERR,"Fix adapt/fep pair style does not exist");
 
312
      if (lmp->suffix_enable) {
 
313
        char psuffix[128];
 
314
        strcpy(psuffix,ad->pstyle);
 
315
        strcat(psuffix,"/");
 
316
        strcat(psuffix,lmp->suffix);
 
317
        pair = force->pair_match(psuffix,1);
 
318
      }
 
319
      if (pair == NULL) pair = force->pair_match(ad->pstyle,1);
 
320
      if (pair == NULL)
 
321
        error->all(FLERR, "Fix adapt/fep pair style does not exist");
223
322
      void *ptr = pair->extract(ad->pparam,ad->pdim);
224
323
      if (ptr == NULL) 
225
324
        error->all(FLERR,"Fix adapt/fep pair style param not supported");
268
367
    }
269
368
  }
270
369
 
271
 
#ifdef ADAPT_DEBUG
272
 
  if (comm->me == 0 && screen) {
273
 
    for (int m = 0; m < nadapt; m++) {
274
 
      Adapt *ad = &adapt[m];
275
 
      if (ad->which == PAIR && ad->pdim == 2) {
276
 
        fprintf(screen, "###ADAPT original %s %s\n", ad->pstyle, ad->pparam);
277
 
        fprintf(screen, "###ADAPT  I  J   old_param\n");
278
 
        for (i = ad->ilo; i <= ad->ihi; i++)
279
 
          for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
280
 
            fprintf(screen, "###ADAPT %2d %2d %9.5f\n", i, j,
281
 
                    ad->array_orig[i][j]);
282
 
      }
283
 
    }
284
 
  }
285
 
#endif
 
370
  // fixes that store initial per-atom values
 
371
  
 
372
  if (id_fix_diam) {
 
373
    int ifix = modify->find_fix(id_fix_diam);
 
374
    if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID");
 
375
    fix_diam = (FixStore *) modify->fix[ifix];
 
376
  }
 
377
  if (id_fix_chg) {
 
378
    int ifix = modify->find_fix(id_fix_chg);
 
379
    if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID");
 
380
    fix_chg = (FixStore *) modify->fix[ifix];
 
381
  }
286
382
 
 
383
  if (strstr(update->integrate_style,"respa"))
 
384
    nlevels_respa = ((Respa *) update->integrate)->nlevels;
287
385
}
288
386
 
289
387
/* ---------------------------------------------------------------------- */
295
393
 
296
394
/* ---------------------------------------------------------------------- */
297
395
 
 
396
void FixAdaptFEP::setup_pre_force_respa(int vflag, int ilevel)
 
397
{
 
398
  if (ilevel < nlevels_respa-1) return;
 
399
  setup_pre_force(vflag);
 
400
}
 
401
 
 
402
/* ---------------------------------------------------------------------- */
 
403
 
298
404
void FixAdaptFEP::pre_force(int vflag)
299
405
{
300
406
  if (nevery == 0) return;
313
419
 
314
420
/* ---------------------------------------------------------------------- */
315
421
 
 
422
void FixAdaptFEP::pre_force_respa(int vflag, int ilevel, int)
 
423
{
 
424
  if (ilevel < nlevels_respa-1) return;
 
425
  pre_force(vflag);
 
426
}
 
427
 
 
428
/* ---------------------------------------------------------------------- */
 
429
 
316
430
void FixAdaptFEP::post_run()
317
431
{
318
432
  if (resetflag) restore_settings();
349
463
          for (i = ad->ilo; i <= ad->ihi; i++)
350
464
            for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
351
465
              ad->array[i][j] = value;
352
 
 
353
 
#ifdef ADAPT_DEBUG
354
 
        if (comm->me == 0 && screen) {
355
 
          fprintf(screen, "###ADAPT change %s %s\n", ad->pstyle, ad->pparam);
356
 
          fprintf(screen, "###ADAPT  I  J   param\n");
357
 
          for (i = ad->ilo; i <= ad->ihi; i++)
358
 
            for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
359
 
              fprintf(screen, "###ADAPT %2d %2d %9.5f\n", i, j, ad->array[i][j]);
360
 
        }
361
 
#endif
362
 
 
363
466
      }
364
467
 
365
468
    // set kspace scale factor
367
470
    } else if (ad->which == KSPACE) {
368
471
      *kspace_scale = value;
369
472
 
 
473
    // set per atom values, also make changes for ghost atoms
 
474
 
370
475
    } else if (ad->which == ATOM) {
371
476
 
372
 
      // set radius from diameter
 
477
      // reset radius from diameter
373
478
      // also scale rmass to new value
374
479
 
375
480
      if (ad->aparam == DIAMETER) {
381
486
        double *radius = atom->radius;
382
487
        double *rmass = atom->rmass;
383
488
        int *mask = atom->mask;
384
 
        int natom = atom->nlocal + atom->nghost;
 
489
        int nlocal = atom->nlocal;
 
490
        int nall = nlocal + atom->nghost;
385
491
 
386
492
        if (mflag == 0) {
387
 
          for (i = 0; i < natom; i++)
 
493
          for (i = 0; i < nall; i++)
388
494
            if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
389
495
              if (mask[i] & groupbit)
390
496
                radius[i] = 0.5*value;
391
497
        } else {
392
 
          for (i = 0; i < natom; i++)
 
498
          for (i = 0; i < nall; i++)
393
499
            if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
394
500
              if (mask[i] & groupbit) {
395
501
                density = rmass[i] / (4.0*MY_PI/3.0 *
403
509
        int *atype = atom->type;
404
510
        double *q = atom->q; 
405
511
        int *mask = atom->mask;
406
 
        int natom = atom->nlocal + atom->nghost;
 
512
        int nlocal = atom->nlocal;
 
513
        int nall = nlocal + atom->nghost;
407
514
 
408
 
        for (i = 0; i < natom; i++)
 
515
        for (i = 0; i < nall; i++)
409
516
          if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
410
517
            if (mask[i] & groupbit) q[i] = value;
411
 
 
412
 
#ifdef ADAPT_DEBUG
413
 
        if (comm->me == 0 && screen) {
414
 
          fprintf(screen, "###ADAPT change charge\n");
415
 
          fprintf(screen, "###ADAPT  atom  I   q\n");
416
 
          for (i = 0; i < atom->nlocal; i++)
417
 
            if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
418
 
              if (mask[i] & groupbit)
419
 
                fprintf(screen, "###ADAPT %5d %2d %9.5f\n", i, atype[i], q[i]);
420
 
        }
421
 
#endif
422
518
      }
423
519
    }
424
520
  }
430
526
  // and also offset and tail corrections
431
527
 
432
528
  if (anypair) force->pair->reinit();
 
529
 
 
530
  // reset KSpace charges if charges have changed
 
531
 
 
532
  if (chgflag && force->kspace) force->kspace->qsum_qsq();
433
533
}
434
534
 
435
535
/* ----------------------------------------------------------------------
452
552
      *kspace_scale = 1.0;
453
553
 
454
554
    } else if (ad->which == ATOM) {
455
 
 
 
555
      if (diamflag) {
 
556
        double density;
 
557
 
 
558
        double *vec = fix_diam->vstore;
 
559
        double *radius = atom->radius;
 
560
        double *rmass = atom->rmass;
 
561
        int *mask = atom->mask;
 
562
        int nlocal = atom->nlocal;
 
563
 
 
564
        for (int i = 0; i < nlocal; i++)
 
565
          if (mask[i] & groupbit) {
 
566
            density = rmass[i] / (4.0*MY_PI/3.0 *
 
567
                                  radius[i]*radius[i]*radius[i]);
 
568
            radius[i] = vec[i];
 
569
            rmass[i] = 4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i] * density;
 
570
          }
 
571
      }
 
572
      if (chgflag) {
 
573
        double *vec = fix_chg->vstore;
 
574
        double *q = atom->q;
 
575
        int *mask = atom->mask;
 
576
        int nlocal = atom->nlocal;
 
577
 
 
578
        for (int i = 0; i < nlocal; i++)
 
579
          if (mask[i] & groupbit) q[i] = vec[i];
 
580
      }
456
581
    }
457
582
  }
458
583
 
459
584
  if (anypair) force->pair->reinit();
 
585
  if (chgflag && force->kspace) force->kspace->qsum_qsq();
 
586
}
 
587
 
 
588
/* ----------------------------------------------------------------------
 
589
   initialize one atom's storage values, called when atom is created
 
590
------------------------------------------------------------------------- */
 
591
 
 
592
void FixAdaptFEP::set_arrays(int i)
 
593
{
 
594
  if (fix_diam) fix_diam->vstore[i] = atom->radius[i];
 
595
  if (fix_chg) fix_chg->vstore[i] = atom->q[i];
460
596
}