← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/USER-FEP/compute_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/compute_fep.cpp
28
28
#include "pair_hybrid.h"
29
29
#include "kspace.h"
30
30
#include "input.h"
 
31
#include "fix.h"
 
32
#include "modify.h"
31
33
#include "variable.h"
32
34
#include "timer.h"
33
35
#include "memory.h"
39
41
enum{PAIR,ATOM};
40
42
enum{CHARGE};
41
43
 
42
 
#undef FEP_DEBUG
43
 
#undef FEP_MAXDEBUG
44
 
 
45
44
/* ---------------------------------------------------------------------- */
46
45
 
47
46
ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
66
65
  int iarg = 4;
67
66
  while (iarg < narg) {
68
67
    if (strcmp(arg[iarg],"pair") == 0) {
69
 
      if (iarg+6 > narg) error->all(FLERR,"Illegal pair attribute in compute fep");
 
68
      if (iarg+6 > narg) error->all(FLERR,
 
69
                                    "Illegal pair attribute in compute fep");
70
70
      npert++;
71
71
      iarg += 6;
72
72
    } else if (strcmp(arg[iarg],"atom") == 0) {
73
 
      if (iarg+4 > narg) error->all(FLERR,"Illegal atom attribute in compute fep");
 
73
      if (iarg+4 > narg) error->all(FLERR,
 
74
                                    "Illegal atom attribute in compute fep");
74
75
      npert++;
75
76
      iarg += 4;
76
77
    } else break;
159
160
 
160
161
  allocate_storage();
161
162
 
 
163
  fixgpu = NULL;
162
164
}
163
165
 
164
166
/* ---------------------------------------------------------------------- */
244
246
                 "compute tail corrections");
245
247
  }
246
248
 
 
249
  // detect if package gpu is present
 
250
 
 
251
  int ifixgpu = modify->find_fix("package_gpu");
 
252
  if (ifixgpu >= 0) fixgpu = modify->fix[ifixgpu];
 
253
 
247
254
  if (comm->me == 0) {
248
255
    if (screen) {
249
256
      fprintf(screen, "FEP settings ...\n");
282
289
{
283
290
  double pe0,pe1;
284
291
 
 
292
  eflag = 1;
 
293
  vflag = 0;
 
294
 
285
295
  invoked_vector = update->ntimestep;
286
296
 
287
297
  if (atom->nmax > nmax) {  // reallocate working arrays if necessary
294
304
 
295
305
  timer->stamp();
296
306
  if (force->pair && force->pair->compute_flag) {
297
 
    force->pair->compute(1,0);
 
307
    force->pair->compute(eflag,vflag);
298
308
    timer->stamp(TIME_PAIR);
299
309
  }
300
 
  if (force->kspace && force->kspace->compute_flag) {
301
 
    force->kspace->compute(1,0);
 
310
  if (chgflag && force->kspace && force->kspace->compute_flag) {
 
311
    force->kspace->compute(eflag,vflag);
302
312
    timer->stamp(TIME_KSPACE);
303
313
  }
 
314
 
 
315
  // accumulate force/energy/virial from /gpu pair styles
 
316
  if (fixgpu) fixgpu->post_force(vflag);
 
317
 
304
318
  pe0 = compute_epair();
305
319
 
306
320
  perturb_params();
307
321
 
308
322
  timer->stamp();
309
323
  if (force->pair && force->pair->compute_flag) {
310
 
    force->pair->compute(1,0);
 
324
    force->pair->compute(eflag,vflag);
311
325
    timer->stamp(TIME_PAIR);
312
326
  }
313
 
  if (force->kspace && force->kspace->compute_flag) {
314
 
    force->kspace->compute(1,0);
 
327
  if (chgflag && force->kspace && force->kspace->compute_flag) {
 
328
    force->kspace->compute(eflag,vflag);
315
329
    timer->stamp(TIME_KSPACE);
316
330
  }
 
331
 
 
332
  // accumulate force/energy/virial from /gpu pair styles
 
333
  // this is required as to empty the answer queue, 
 
334
  // otherwise the force compute on the GPU in the next step would be incorrect
 
335
  if (fixgpu) fixgpu->post_force(vflag);
 
336
    
317
337
  pe1 = compute_epair();
318
338
 
319
339
  restore_qfev();   // restore charge, force, energy, virial array values
324
344
  vector[2] = domain->xprd * domain->yprd * domain->zprd;
325
345
  if (volumeflag)
326
346
    vector[1] *= vector[2];
327
 
 
328
 
#ifdef FEP_DEBUG
329
 
  if (comm->me == 0 && screen)
330
 
    fprintf(screen, "FEP U0 = %f  U1 = %f  DU = %f  exp(-DU/kT) = %f\n",
331
 
            pe0,pe1,vector[0],vector[1]);
332
 
#endif
333
347
}
334
348
 
335
349
 
351
365
    eng_pair += force->pair->etail / volume;
352
366
  }
353
367
 
354
 
  if (force->kspace) eng_pair += force->kspace->energy;
 
368
  if (chgflag && force->kspace) eng_pair += force->kspace->energy;
355
369
 
356
370
  return eng_pair;
357
371
}
374
388
      for (i = pert->ilo; i <= pert->ihi; i++)
375
389
        for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
376
390
          pert->array[i][j] = pert->array_orig[i][j] + delta;
377
 
      
378
 
#ifdef FEP_MAXDEBUG
379
 
      if (comm->me == 0 && screen) {
380
 
        fprintf(screen, "###FEP change %s %s, delta = %f\n",
381
 
                pert->pstyle, pert->pparam, delta);
382
 
        fprintf(screen, "###FEP  I  J   old_param new_param\n");
383
 
        for (i = pert->ilo; i <= pert->ihi; i++)
384
 
          for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
385
 
            fprintf(screen, "###FEP %2d %2d %9.5f %9.5f\n", i, j, 
386
 
                    pert->array_orig[i][j], pert->array[i][j]);
387
 
      }
388
 
#endif
389
391
 
390
392
    } else if (pert->which == ATOM) {
391
393
 
400
402
            if (mask[i] & groupbit)
401
403
              q[i] += delta; 
402
404
 
403
 
#ifdef FEP_MAXDEBUG
404
 
        if (comm->me == 0 && screen) {
405
 
          fprintf(screen, "###FEP change charge, delta = %f\n", delta);
406
 
          fprintf(screen, "###FEP  atom  I   old_q     new_q\n");
407
 
          for (i = 0; i < atom->nlocal; i++)
408
 
            if (atype[i] >= pert->ilo && atype[i] <= pert->ihi)
409
 
              if (mask[i] & groupbit)
410
 
                fprintf(screen, "###FEP %5d %2d %9.5f %9.5f\n", i, atype[i],
411
 
                        q_orig[i], q[i]);
412
 
        }
413
 
#endif
414
 
 
415
405
      }
416
406
    }
417
407
  }
421
411
  // and also offset and tail corrections
422
412
 
423
413
  if (pairflag) force->pair->reinit();
 
414
 
 
415
  // reset KSpace charges if charges have changed
 
416
 
 
417
  if (chgflag && force->kspace) force->kspace->qsum_qsq();
424
418
}
425
419
 
426
420
 
457
451
      for (i = pert->ilo; i <= pert->ihi; i++)
458
452
        for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
459
453
          pert->array[i][j] = pert->array_orig[i][j];
460
 
 
461
 
#ifdef FEP_MAXDEBUG
462
 
      if (comm->me == 0 && screen) {
463
 
        fprintf(screen, "###FEP restore %s %s\n", pert->pstyle, pert->pparam);
464
 
        fprintf(screen, "###FEP  I  J   param\n");
465
 
        for (i = pert->ilo; i <= pert->ihi; i++)
466
 
          for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
467
 
            fprintf(screen, "###FEP %2d %2d %9.5f\n", i, j, pert->array[i][j]);
468
 
      }
469
 
#endif
470
 
    }
471
 
#ifdef FEP_MAXDEBUG
472
 
    if (pert->which == ATOM) {
473
 
      if (comm->me == 0 && screen) {
474
 
        int *atype = atom->type;
475
 
        double *q = atom->q; 
476
 
        int *mask = atom->mask;
477
 
        int natom = atom->nlocal;
478
 
        fprintf(screen, "###FEP restore charge\n");
479
 
        fprintf(screen, "###FEP  atom  I   q\n");
480
 
        for (i = 0; i < natom; i++)
481
 
          if (atype[i] >= pert->ilo && atype[i] <= pert->ihi)
482
 
            if (mask[i] & groupbit)
483
 
              fprintf(screen, "###FEP %5d %2d %9.5f\n", i, atype[i], q[i]);
484
 
      }
485
 
    }
486
 
#endif
 
454
    }
487
455
  }
488
456
 
489
 
  // re-initialize pair styles if any PAIR settings were changed
490
 
  // this resets other coeffs that may depend on changed values,
491
 
  // and also offset and tail corrections
492
 
 
493
457
  if (pairflag) force->pair->reinit();
 
458
 
 
459
  // reset KSpace charges if charges have changed
 
460
 
 
461
  if (chgflag && force->kspace) force->kspace->qsum_qsq();
494
462
}
495
463
 
496
464
 
501
469
void ComputeFEP::allocate_storage()
502
470
{
503
471
  nmax = atom->nmax;
504
 
  if (chgflag)
505
 
    memory->create(q_orig,nmax,"fep:q_orig");
506
472
  memory->create(f_orig,nmax,3,"fep:f_orig");
507
473
  memory->create(peatom_orig,nmax,"fep:peatom_orig");
508
474
  memory->create(pvatom_orig,nmax,6,"fep:pvatom_orig");
509
 
  if (force->kspace) {
510
 
    memory->create(keatom_orig,nmax,"fep:keatom_orig");
511
 
    memory->create(kvatom_orig,nmax,6,"fep:kvatom_orig");
 
475
  if (chgflag) {
 
476
    memory->create(q_orig,nmax,"fep:q_orig");
 
477
    if (force->kspace) {
 
478
      memory->create(keatom_orig,nmax,"fep:keatom_orig");
 
479
      memory->create(kvatom_orig,nmax,6,"fep:kvatom_orig");
 
480
    }
512
481
  }
513
482
}
514
483
 
516
485
 
517
486
void ComputeFEP::deallocate_storage()
518
487
{
519
 
  if (chgflag)
520
 
    memory->destroy(q_orig);
521
488
  memory->destroy(f_orig);
522
489
  memory->destroy(peatom_orig);
523
490
  memory->destroy(pvatom_orig);
524
 
  if (force->kspace) {
525
 
    memory->destroy(keatom_orig);
526
 
    memory->destroy(kvatom_orig);
 
491
  if (chgflag) {
 
492
    memory->destroy(q_orig);
 
493
    if (force->kspace) {
 
494
      memory->destroy(keatom_orig);
 
495
      memory->destroy(kvatom_orig);
 
496
    }
527
497
  }
528
498
}
529
499
 
541
511
  if (force->newton || force->kspace->tip4pflag)
542
512
      natom += atom->nghost;
543
513
 
544
 
  if (chgflag) {
545
 
    double *q = atom->q; 
546
 
    for (i = 0; i < nall; i++)
547
 
      q_orig[i] = q[i];
548
 
  }
549
 
 
550
514
  double **f = atom->f;
551
515
  for (i = 0; i < natom; i++) {
552
516
    f_orig[i][0] = f[i][0]; 
581
545
    }
582
546
  }
583
547
 
584
 
  if (force->kspace) {
585
 
    energy_orig = force->kspace->energy;
586
 
    kvirial_orig[0] = force->kspace->virial[0];
587
 
    kvirial_orig[1] = force->kspace->virial[1];
588
 
    kvirial_orig[2] = force->kspace->virial[2];
589
 
    kvirial_orig[3] = force->kspace->virial[3];
590
 
    kvirial_orig[4] = force->kspace->virial[4];
591
 
    kvirial_orig[5] = force->kspace->virial[5];
592
 
    
593
 
    if (update->eflag_atom) {
594
 
      double *keatom = force->kspace->eatom;
595
 
      for (i = 0; i < natom; i++)
596
 
        keatom_orig[i] = keatom[i];
597
 
    }
598
 
    if (update->vflag_atom) {
599
 
      double **kvatom = force->kspace->vatom;
600
 
      for (i = 0; i < natom; i++) {
601
 
        kvatom_orig[i][0] = kvatom[i][0];
602
 
        kvatom_orig[i][1] = kvatom[i][1];
603
 
        kvatom_orig[i][2] = kvatom[i][2];
604
 
        kvatom_orig[i][3] = kvatom[i][3];
605
 
        kvatom_orig[i][4] = kvatom[i][4];
606
 
        kvatom_orig[i][5] = kvatom[i][5];
 
548
  if (chgflag) {
 
549
    double *q = atom->q; 
 
550
    for (i = 0; i < nall; i++)
 
551
      q_orig[i] = q[i];
 
552
 
 
553
    if (force->kspace) {
 
554
      energy_orig = force->kspace->energy;
 
555
      kvirial_orig[0] = force->kspace->virial[0];
 
556
      kvirial_orig[1] = force->kspace->virial[1];
 
557
      kvirial_orig[2] = force->kspace->virial[2];
 
558
      kvirial_orig[3] = force->kspace->virial[3];
 
559
      kvirial_orig[4] = force->kspace->virial[4];
 
560
      kvirial_orig[5] = force->kspace->virial[5];
 
561
      
 
562
      if (update->eflag_atom) {
 
563
        double *keatom = force->kspace->eatom;
 
564
        for (i = 0; i < natom; i++)
 
565
          keatom_orig[i] = keatom[i];
 
566
      }
 
567
      if (update->vflag_atom) {
 
568
        double **kvatom = force->kspace->vatom;
 
569
        for (i = 0; i < natom; i++) {
 
570
          kvatom_orig[i][0] = kvatom[i][0];
 
571
          kvatom_orig[i][1] = kvatom[i][1];
 
572
          kvatom_orig[i][2] = kvatom[i][2];
 
573
          kvatom_orig[i][3] = kvatom[i][3];
 
574
          kvatom_orig[i][4] = kvatom[i][4];
 
575
          kvatom_orig[i][5] = kvatom[i][5];
 
576
        }
607
577
      }
608
578
    }
609
579
  }
620
590
  if (force->newton || force->kspace->tip4pflag)
621
591
      natom += atom->nghost;
622
592
 
623
 
  if (chgflag) {
624
 
    double *q = atom->q; 
625
 
    for (i = 0; i < nall; i++)
626
 
      q[i] = q_orig[i];
627
 
  }
628
 
 
629
593
  double **f = atom->f;
630
594
  for (i = 0; i < natom; i++) {
631
595
    f[i][0] = f_orig[i][0]; 
660
624
    }
661
625
  }
662
626
 
663
 
  if (force->kspace) {
664
 
    force->kspace->energy = energy_orig;
665
 
    force->kspace->virial[0] = kvirial_orig[0];
666
 
    force->kspace->virial[1] = kvirial_orig[1];
667
 
    force->kspace->virial[2] = kvirial_orig[2];
668
 
    force->kspace->virial[3] = kvirial_orig[3];
669
 
    force->kspace->virial[4] = kvirial_orig[4];
670
 
    force->kspace->virial[5] = kvirial_orig[5];
 
627
  if (chgflag) {
 
628
    double *q = atom->q; 
 
629
    for (i = 0; i < nall; i++)
 
630
      q[i] = q_orig[i];
 
631
 
 
632
    if (force->kspace) {
 
633
      force->kspace->energy = energy_orig;
 
634
      force->kspace->virial[0] = kvirial_orig[0];
 
635
      force->kspace->virial[1] = kvirial_orig[1];
 
636
      force->kspace->virial[2] = kvirial_orig[2];
 
637
      force->kspace->virial[3] = kvirial_orig[3];
 
638
      force->kspace->virial[4] = kvirial_orig[4];
 
639
      force->kspace->virial[5] = kvirial_orig[5];
671
640
    
672
 
    if (update->eflag_atom) {
673
 
      double *keatom = force->kspace->eatom;
674
 
      for (i = 0; i < natom; i++)
675
 
        keatom[i] = keatom_orig[i];
676
 
    }
677
 
    if (update->vflag_atom) {
678
 
      double **kvatom = force->kspace->vatom;
679
 
      for (i = 0; i < natom; i++) {
680
 
        kvatom[i][0] = kvatom_orig[i][0];
681
 
        kvatom[i][1] = kvatom_orig[i][1];
682
 
        kvatom[i][2] = kvatom_orig[i][2];
683
 
        kvatom[i][3] = kvatom_orig[i][3];
684
 
        kvatom[i][4] = kvatom_orig[i][4];
685
 
        kvatom[i][5] = kvatom_orig[i][5];
 
641
      if (update->eflag_atom) {
 
642
        double *keatom = force->kspace->eatom;
 
643
        for (i = 0; i < natom; i++)
 
644
          keatom[i] = keatom_orig[i];
 
645
      }
 
646
      if (update->vflag_atom) {
 
647
        double **kvatom = force->kspace->vatom;
 
648
        for (i = 0; i < natom; i++) {
 
649
          kvatom[i][0] = kvatom_orig[i][0];
 
650
          kvatom[i][1] = kvatom_orig[i][1];
 
651
          kvatom[i][2] = kvatom_orig[i][2];
 
652
          kvatom[i][3] = kvatom_orig[i][3];
 
653
          kvatom[i][4] = kvatom_orig[i][4];
 
654
          kvatom[i][5] = kvatom_orig[i][5];
 
655
        }
686
656
      }
687
657
    }
688
658
  }