← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/RIGID/fix_rigid.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/RIGID/fix_rigid.cpp
75
75
 
76
76
  extended = orientflag = dorientflag = 0;
77
77
  body = NULL;
 
78
  xcmimage = NULL;
78
79
  displace = NULL;
79
80
  eflags = NULL;
80
81
  orient = NULL;
336
337
 
337
338
    } else if (strcmp(arg[iarg],"langevin") == 0) {
338
339
      if (iarg+5 > narg) error->all(FLERR,"Illegal fix rigid command");
339
 
      if (strcmp(style,"rigid") != 0 && strcmp(style,"rigid/nve") != 0)
 
340
      if (strcmp(style,"rigid") != 0 && strcmp(style,"rigid/nve") != 0 && 
 
341
          strcmp(style,"rigid/omp") != 0 && strcmp(style,"rigid/nve/omp") != 0)
340
342
        error->all(FLERR,"Illegal fix rigid command");
341
343
      langflag = 1;
342
344
      t_start = force->numeric(FLERR,arg[iarg+1]);
350
352
 
351
353
    } else if (strcmp(arg[iarg],"temp") == 0) {
352
354
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
353
 
      if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0)
 
355
      if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0 &&
 
356
          strcmp(style,"rigid/nvt/omp") != 0 && 
 
357
          strcmp(style,"rigid/npt/omp") != 0)
354
358
        error->all(FLERR,"Illegal fix rigid command");
355
359
      tstat_flag = 1;
356
360
      t_start = force->numeric(FLERR,arg[iarg+1]);
360
364
 
361
365
    } else if (strcmp(arg[iarg],"iso") == 0) {
362
366
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
363
 
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0)
364
 
              error->all(FLERR,"Illegal fix rigid command");
 
367
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
 
368
          strcmp(style,"rigid/npt/omp") != 0 && 
 
369
          strcmp(style,"rigid/nph/omp") != 0)
 
370
        error->all(FLERR,"Illegal fix rigid command");
365
371
      pcouple = XYZ;
366
372
      p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
367
373
      p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
376
382
 
377
383
    } else if (strcmp(arg[iarg],"aniso") == 0) {
378
384
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
379
 
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0)
380
 
              error->all(FLERR,"Illegal fix rigid command");
 
385
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
 
386
          strcmp(style,"rigid/npt/omp") != 0 && 
 
387
          strcmp(style,"rigid/nph/omp") != 0)
 
388
        error->all(FLERR,"Illegal fix rigid command");
381
389
      p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
382
390
      p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
383
391
      p_period[0] = p_period[1] = p_period[2] = 
391
399
 
392
400
    } else if (strcmp(arg[iarg],"x") == 0) {
393
401
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
 
402
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
 
403
          strcmp(style,"rigid/npt/omp") != 0 && 
 
404
          strcmp(style,"rigid/nph/omp") != 0)
 
405
        error->all(FLERR,"Illegal fix rigid command");
394
406
      p_start[0] = force->numeric(FLERR,arg[iarg+1]);
395
407
      p_stop[0] = force->numeric(FLERR,arg[iarg+2]);
396
408
      p_period[0] = force->numeric(FLERR,arg[iarg+3]);
399
411
 
400
412
    } else if (strcmp(arg[iarg],"y") == 0) {
401
413
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
 
414
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
 
415
          strcmp(style,"rigid/npt/omp") != 0 && 
 
416
          strcmp(style,"rigid/nph/omp") != 0)
 
417
        error->all(FLERR,"Illegal fix rigid command");
402
418
      p_start[1] = force->numeric(FLERR,arg[iarg+1]);
403
419
      p_stop[1] = force->numeric(FLERR,arg[iarg+2]);
404
420
      p_period[1] = force->numeric(FLERR,arg[iarg+3]);
407
423
 
408
424
    } else if (strcmp(arg[iarg],"z") == 0) {
409
425
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
 
426
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
 
427
          strcmp(style,"rigid/npt/omp") != 0 && 
 
428
          strcmp(style,"rigid/nph/omp") != 0)
 
429
        error->all(FLERR,"Illegal fix rigid command");
410
430
      p_start[2] = force->numeric(FLERR,arg[iarg+1]);
411
431
      p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
412
432
      p_period[2] = force->numeric(FLERR,arg[iarg+3]);
425
445
 
426
446
    } else if (strcmp(arg[iarg],"dilate") == 0) {
427
447
      if (iarg+2 > narg) 
428
 
        error->all(FLERR,"Illegal fix rigid nvt/npt/nph command");
 
448
        error->all(FLERR,"Illegal fix rigid npt/nph command");
429
449
      if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
430
450
      else {
431
451
        allremap = 0;
442
462
 
443
463
    } else if (strcmp(arg[iarg],"tparam") == 0) {
444
464
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
445
 
      if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0)
 
465
      if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0 && 
 
466
          strcmp(style,"rigid/nvt/omp") != 0 && 
 
467
          strcmp(style,"rigid/npt/omp") != 0)
446
468
        error->all(FLERR,"Illegal fix rigid command");
447
469
      t_chain = force->inumeric(FLERR,arg[iarg+1]);
448
470
      t_iter = force->inumeric(FLERR,arg[iarg+2]);
451
473
 
452
474
    } else if (strcmp(arg[iarg],"pchain") == 0) {
453
475
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
454
 
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0)
 
476
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 && 
 
477
          strcmp(style,"rigid/npt/omp") != 0 && 
 
478
          strcmp(style,"rigid/nph/omp") != 0)
455
479
        error->all(FLERR,"Illegal fix rigid command");
456
480
      p_chain = force->inumeric(FLERR,arg[iarg+1]);
457
481
      iarg += 2;
523
547
  MINUSPI = -MY_PI;
524
548
  TWOPI = 2.0*MY_PI;
525
549
 
 
550
  // if infile set, only setup rigid bodies once, using info from file
 
551
  // this means users cannot change atom properties like mass between runs
 
552
  // if they do, bodies will not reflect the changes
 
553
 
 
554
  staticflag = 0;
 
555
  if (infile) setup_bodies_static();
 
556
 
526
557
  // print statistics
527
558
 
528
559
  int nsum = 0;
532
563
    if (screen) fprintf(screen,"%d rigid bodies with %d atoms\n",nbody,nsum);
533
564
    if (logfile) fprintf(logfile,"%d rigid bodies with %d atoms\n",nbody,nsum);
534
565
  }
535
 
 
536
 
  // firstflag = 1 triggers one-time initialization of rigid body attributes
537
 
 
538
 
  firstflag = 1;
539
566
}
540
567
 
541
568
/* ---------------------------------------------------------------------- */
551
578
  memory->destroy(mol2body);
552
579
  memory->destroy(body2mol);
553
580
 
554
 
  // delete locally stored arrays
 
581
  // delete locally stored per-atom arrays
555
582
 
556
583
  memory->destroy(body);
 
584
  memory->destroy(xcmimage);
557
585
  memory->destroy(displace);
558
586
  memory->destroy(eflags);
559
587
  memory->destroy(orient);
640
668
  if (strstr(update->integrate_style,"respa"))
641
669
    step_respa = ((Respa *) update->integrate)->step;
642
670
 
643
 
  // one-time initialization of rigid body attributes
644
 
  // setup_bodies_static = extended flags, masstotal, COM, inertia tensor
645
 
  // setup_bodies_dynamic = vcm and angmom
 
671
  // setup rigid bodies, using current atom info
 
672
  // allows resetting of atom properties like mass between runs
 
673
  // only do this if not using an infile, else was called in constructor
646
674
 
647
 
  if (firstflag) {
648
 
    firstflag = 0;
649
 
    setup_bodies_static();
650
 
    setup_bodies_dynamic();
651
 
  }
 
675
  if (!infile) setup_bodies_static();
652
676
 
653
677
  // temperature scale factor
654
678
 
662
686
}
663
687
 
664
688
/* ----------------------------------------------------------------------
 
689
   invoke pre_neighbor() to insure body xcmimage flags are reset
 
690
     needed if Verlet::setup::pbc() has remapped/migrated atoms for 2nd run
 
691
------------------------------------------------------------------------- */
 
692
 
 
693
void FixRigid::setup_pre_neighbor()
 
694
{
 
695
  pre_neighbor();
 
696
}
 
697
 
 
698
/* ----------------------------------------------------------------------
665
699
   compute initial fcm and torque on bodies, also initial virial
666
700
   reset all particle velocities to be consistent with vcm and omega
667
701
------------------------------------------------------------------------- */
670
704
{
671
705
  int i,n,ibody;
672
706
 
 
707
  // setup_bodies_dynamic sets vcm and angmom
 
708
  // so angmom_to_omega() and set_v() below will set per-atom vels correctly
 
709
  // re-calling it every run allows reset of body/atom velocities between runs
 
710
 
 
711
  setup_bodies_dynamic();
 
712
 
673
713
  // fcm = force on center-of-mass of each rigid body
674
714
 
675
715
  double **f = atom->f;
696
736
 
697
737
  // torque = torque on each rigid body
698
738
 
699
 
  imageint *image = atom->image;
700
739
  double **x = atom->x;
701
740
 
702
741
  double dx,dy,dz;
709
748
    if (body[i] < 0) continue;
710
749
    ibody = body[i];
711
750
 
712
 
    domain->unmap(x[i],image[i],unwrap);
 
751
    domain->unmap(x[i],xcmimage[i],unwrap);
713
752
    dx = unwrap[0] - xcm[ibody][0];
714
753
    dy = unwrap[1] - xcm[ibody][1];
715
754
    dz = unwrap[2] - xcm[ibody][2];
760
799
  for (ibody = 0; ibody < nbody; ibody++)
761
800
    MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
762
801
                               ez_space[ibody],inertia[ibody],omega[ibody]);
 
802
 
763
803
  set_v();
764
804
 
765
805
  // guesstimate virial as 2x the set_v contribution
877
917
 
878
918
  // sum over atoms to get force and torque on rigid body
879
919
 
880
 
  imageint *image = atom->image;
881
920
  double **x = atom->x;
882
921
  double **f = atom->f;
883
922
  int nlocal = atom->nlocal;
896
935
    sum[ibody][1] += f[i][1];
897
936
    sum[ibody][2] += f[i][2];
898
937
 
899
 
    domain->unmap(x[i],image[i],unwrap);
 
938
    domain->unmap(x[i],xcmimage[i],unwrap);
900
939
    dx = unwrap[0] - xcm[ibody][0];
901
940
    dy = unwrap[1] - xcm[ibody][1];
902
941
    dz = unwrap[2] - xcm[ibody][2];
960
999
  set_v();
961
1000
}
962
1001
 
963
 
/* ----------------------------------------------------------------------
964
 
   apply evolution operators to quat, quat momentum
965
 
   see Miller paper cited in fix rigid/nvt and fix rigid/npt
966
 
------------------------------------------------------------------------- */
967
 
 
968
 
void FixRigid::no_squish_rotate(int k, double *p, double *q,
969
 
                                double *inertia, double dt) const
970
 
{
971
 
  double phi,c_phi,s_phi,kp[4],kq[4];
972
 
 
973
 
  // apply permuation operator on p and q, get kp and kq
974
 
 
975
 
  if (k == 1) {
976
 
    kq[0] = -q[1];  kp[0] = -p[1];
977
 
    kq[1] =  q[0];  kp[1] =  p[0];
978
 
    kq[2] =  q[3];  kp[2] =  p[3];
979
 
    kq[3] = -q[2];  kp[3] = -p[2];
980
 
  } else if (k == 2) {
981
 
    kq[0] = -q[2];  kp[0] = -p[2];
982
 
    kq[1] = -q[3];  kp[1] = -p[3];
983
 
    kq[2] =  q[0];  kp[2] =  p[0];
984
 
    kq[3] =  q[1];  kp[3] =  p[1];
985
 
  } else if (k == 3) {
986
 
    kq[0] = -q[3];  kp[0] = -p[3];
987
 
    kq[1] =  q[2];  kp[1] =  p[2];
988
 
    kq[2] = -q[1];  kp[2] = -p[1];
989
 
    kq[3] =  q[0];  kp[3] =  p[0];
990
 
  }
991
 
 
992
 
  // obtain phi, cosines and sines
993
 
 
994
 
  phi = p[0]*kq[0] + p[1]*kq[1] + p[2]*kq[2] + p[3]*kq[3];
995
 
  if (fabs(inertia[k-1]) < 1e-6) phi *= 0.0;
996
 
  else phi /= 4.0 * inertia[k-1];
997
 
  c_phi = cos(dt * phi);
998
 
  s_phi = sin(dt * phi);
999
 
 
1000
 
  // advance p and q
1001
 
 
1002
 
  p[0] = c_phi*p[0] + s_phi*kp[0];
1003
 
  p[1] = c_phi*p[1] + s_phi*kp[1];
1004
 
  p[2] = c_phi*p[2] + s_phi*kp[2];
1005
 
  p[3] = c_phi*p[3] + s_phi*kp[3];
1006
 
 
1007
 
  q[0] = c_phi*q[0] + s_phi*kq[0];
1008
 
  q[1] = c_phi*q[1] + s_phi*kq[1];
1009
 
  q[2] = c_phi*q[2] + s_phi*kq[2];
1010
 
  q[3] = c_phi*q[3] + s_phi*kq[3];
1011
 
}
1012
 
 
1013
1002
/* ---------------------------------------------------------------------- */
1014
1003
 
1015
1004
void FixRigid::initial_integrate_respa(int vflag, int ilevel, int iloop)
1035
1024
   done during pre_neighbor so will be after call to pbc()
1036
1025
     and after fix_deform::pre_exchange() may have flipped box
1037
1026
   use domain->remap() in case xcm is far away from box
1038
 
     due to 1st definition of rigid body or due to box flip
1039
 
   if don't do this, then atoms of a body which drifts far away
1040
 
     from a triclinic box will be remapped back into box
1041
 
     with huge displacements when the box tilt changes via set_x()
1042
 
   adjust image flag of body and image flags of all atoms in body
 
1027
     due to first-time definition of rigid body in setup_bodies_static()
 
1028
     or due to box flip
 
1029
   also adjust imagebody = rigid body image flags, due to xcm remap
 
1030
   also reset body xcmimage flags of all atoms in bodies
 
1031
   xcmimage flags are relative to xcm so that body can be unwrapped
 
1032
   if don't do this, would need xcm to move with true image flags
 
1033
     then a body could end up very far away from box
 
1034
     set_xv() will then compute huge displacements every step to 
 
1035
       reset coords of all body atoms to be back inside the box,
 
1036
       ditto for triclinic box flip, which causes numeric problems
1043
1037
------------------------------------------------------------------------- */
1044
1038
 
1045
1039
void FixRigid::pre_neighbor()
1046
1040
{
1047
 
  imageint original,oldimage,newimage;
1048
 
 
1049
 
  for (int ibody = 0; ibody < nbody; ibody++) {
1050
 
    original = imagebody[ibody];
 
1041
  for (int ibody = 0; ibody < nbody; ibody++)
1051
1042
    domain->remap(xcm[ibody],imagebody[ibody]);
1052
 
 
1053
 
    if (original == imagebody[ibody]) remapflag[ibody][3] = 0;
1054
 
    else {
1055
 
      oldimage = original & IMGMASK;
1056
 
      newimage = imagebody[ibody] & IMGMASK;
1057
 
      remapflag[ibody][0] = newimage - oldimage;
1058
 
      oldimage = (original >> IMGBITS) & IMGMASK;
1059
 
      newimage = (imagebody[ibody] >> IMGBITS) & IMGMASK;
1060
 
      remapflag[ibody][1] = newimage - oldimage;
1061
 
      oldimage = original >> IMG2BITS;
1062
 
      newimage = imagebody[ibody] >> IMG2BITS;
1063
 
      remapflag[ibody][2] = newimage - oldimage;
1064
 
      remapflag[ibody][3] = 1;
1065
 
    }
1066
 
  }
1067
 
 
1068
 
  // adjust image flags of any atom in a rigid body whose xcm was remapped
1069
 
  // subtracting remapflag = new-old keeps ix,iy,iz near 0
1070
 
  //   so body is always in central simulation box
 
1043
  image_shift();
 
1044
}
 
1045
 
 
1046
/* ----------------------------------------------------------------------
 
1047
   reset body xcmimage flags of atoms in bodies
 
1048
   xcmimage flags are relative to xcm so that body can be unwrapped
 
1049
   xcmimage = true image flag - imagebody flag
 
1050
------------------------------------------------------------------------- */
 
1051
 
 
1052
void FixRigid::image_shift()
 
1053
{
 
1054
  int ibody;
 
1055
  imageint tdim,bdim,xdim[3];
1071
1056
 
1072
1057
  imageint *image = atom->image;
1073
1058
  int nlocal = atom->nlocal;
1074
1059
 
1075
 
  int ibody;
1076
 
  imageint idim,otherdims;
1077
 
 
1078
1060
  for (int i = 0; i < nlocal; i++) {
1079
 
    if (body[i] == -1) continue;
1080
 
    if (remapflag[body[i]][3] == 0) continue;
 
1061
    if (body[i] < 0) continue;
1081
1062
    ibody = body[i];
1082
1063
 
1083
 
    if (remapflag[ibody][0]) {
1084
 
      idim = image[i] & IMGMASK;
1085
 
      otherdims = image[i] ^ idim;
1086
 
      idim -= remapflag[ibody][0];
1087
 
      idim &= IMGMASK;
1088
 
      image[i] = otherdims | idim;
1089
 
    }
1090
 
    if (remapflag[ibody][1]) {
1091
 
      idim = (image[i] >> IMGBITS) & IMGMASK;
1092
 
      otherdims = image[i] ^ (idim << IMGBITS);
1093
 
      idim -= remapflag[ibody][1];
1094
 
      idim &= IMGMASK;
1095
 
      image[i] = otherdims | (idim << IMGBITS);
1096
 
    }
1097
 
    if (remapflag[ibody][2]) {
1098
 
      idim = image[i] >> IMG2BITS;
1099
 
      otherdims = image[i] ^ (idim << IMG2BITS);
1100
 
      idim -= remapflag[ibody][2];
1101
 
      idim &= IMGMASK;
1102
 
      image[i] = otherdims | (idim << IMG2BITS);
1103
 
    }
 
1064
    tdim = image[i] & IMGMASK;
 
1065
    bdim = imagebody[ibody] & IMGMASK;
 
1066
    xdim[0] = IMGMAX + tdim - bdim;
 
1067
    tdim = (image[i] >> IMGBITS) & IMGMASK;
 
1068
    bdim = (imagebody[ibody] >> IMGBITS) & IMGMASK;
 
1069
    xdim[1] = IMGMAX + tdim - bdim;
 
1070
    tdim = image[i] >> IMG2BITS;
 
1071
    bdim = imagebody[ibody] >> IMG2BITS;
 
1072
    xdim[2] = IMGMAX + tdim - bdim;
 
1073
 
 
1074
    xcmimage[i] = (xdim[2] << IMG2BITS) | (xdim[1] << IMGBITS) | xdim[0];
1104
1075
  }
1105
1076
}
1106
1077
 
1113
1084
{
1114
1085
  // cannot count DOF correctly unless setup_bodies_static() has been called
1115
1086
 
1116
 
  if (firstflag) {
 
1087
  if (!staticflag) {
1117
1088
    if (comm->me == 0) 
1118
1089
      error->warning(FLERR,"Cannot count rigid body degrees-of-freedom "
1119
 
                     "before bodies are fully initialized");
 
1090
                     "before bodies are initialized");
1120
1091
    return 0;
1121
1092
  }
1122
1093
 
1136
1107
 
1137
1108
  for (int i = 0; i < nlocal; i++)
1138
1109
    if (body[i] >= 0 && mask[i] & tgroupbit) {
1139
 
      if (extended && eflags[i]) mcount[body[i]]++;
 
1110
      // do not count point particles or point dipoles as extended particles
 
1111
      // a spheroid dipole will be counted as extended
 
1112
      if (extended && (eflags[i] & ~(POINT | DIPOLE))) mcount[body[i]]++;
1140
1113
      else ncount[body[i]]++;
1141
1114
    }
1142
1115
 
1219
1192
  double xy,xz,yz;
1220
1193
  double ione[3],exone[3],eyone[3],ezone[3],vr[6],p[3][3];
1221
1194
 
1222
 
  imageint *image = atom->image;
1223
1195
  double **x = atom->x;
1224
1196
  double **v = atom->v;
1225
1197
  double **f = atom->f;
1244
1216
    if (body[i] < 0) continue;
1245
1217
    ibody = body[i];
1246
1218
 
1247
 
    xbox = (image[i] & IMGMASK) - IMGMAX;
1248
 
    ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
1249
 
    zbox = (image[i] >> IMG2BITS) - IMGMAX;
 
1219
    xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
 
1220
    ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
 
1221
    zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
1250
1222
 
1251
1223
    // save old positions and velocities for virial
1252
1224
 
1401
1373
  double *rmass = atom->rmass;
1402
1374
  double *mass = atom->mass;
1403
1375
  int *type = atom->type;
1404
 
  imageint *image = atom->image;
1405
1376
  int nlocal = atom->nlocal;
1406
1377
 
1407
1378
  double xprd = domain->xprd;
1450
1421
      fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
1451
1422
      fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
1452
1423
 
1453
 
      xbox = (image[i] & IMGMASK) - IMGMAX;
1454
 
      ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
1455
 
      zbox = (image[i] >> IMG2BITS) - IMGMAX;
 
1424
      xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
 
1425
      ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
 
1426
      zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
1456
1427
 
1457
1428
      if (triclinic == 0) {
1458
1429
        x0 = x[i][0] + xbox*xprd;
1522
1493
}
1523
1494
 
1524
1495
/* ----------------------------------------------------------------------
1525
 
   one-time initialization of static rigid body attributes
1526
 
   extended flags, masstotal, center-of-mass
1527
 
   Cartesian and diagonalized inertia tensor
1528
 
   read per-body attributes from infile if specified
 
1496
   initialization of static rigid body attributes
 
1497
   called from init() so body/atom properties can be changed between runs
 
1498
   unless reading from infile, in which case called once from constructor
 
1499
   sets extended flags, masstotal, center-of-mass
 
1500
   sets Cartesian and diagonalized inertia tensor
 
1501
   sets body image flags, but only on first call
1529
1502
------------------------------------------------------------------------- */
1530
1503
 
1531
1504
void FixRigid::setup_bodies_static()
1611
1584
    }
1612
1585
  }
1613
1586
 
 
1587
  // first-time setting of body xcmimage flags = true image flags
 
1588
 
 
1589
  if (!staticflag) {
 
1590
    imageint *image = atom->image;
 
1591
    for (i = 0; i < nlocal; i++)
 
1592
      if (body[i] >= 0) xcmimage[i] = image[i];
 
1593
      else xcmimage[i] = 0;
 
1594
  }
 
1595
 
1614
1596
  // compute masstotal & center-of-mass of each rigid body
1615
1597
  // error if image flag is not 0 in a non-periodic dim
1616
1598
 
1617
1599
  double **x = atom->x;
1618
 
  imageint *image = atom->image;
1619
1600
 
1620
1601
  int *periodicity = domain->periodicity;
1621
1602
  double xprd = domain->xprd;
1634
1615
    if (body[i] < 0) continue;
1635
1616
    ibody = body[i];
1636
1617
 
1637
 
    xbox = (image[i] & IMGMASK) - IMGMAX;
1638
 
    ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
1639
 
    zbox = (image[i] >> IMG2BITS) - IMGMAX;
 
1618
    xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
 
1619
    ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
 
1620
    zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
1640
1621
    if (rmass) massone = rmass[i];
1641
1622
    else massone = mass[type[i]];
1642
1623
 
1680
1661
    readfile(0,masstotal,xcm,inbody);
1681
1662
  }
1682
1663
 
1683
 
  // set image flags for each rigid body to default values
1684
 
  // then remap the xcm of each body back into simulation box if needed
 
1664
  // one-time set of rigid body image flags to default values
 
1665
  //   staticflag insures this is only done once, not on successive runs
 
1666
  // then remap the xcm of each body back into simulation box
 
1667
  //   and reset body xcmimage flags via pre_neighbor()
1685
1668
 
1686
 
  for (ibody = 0; ibody < nbody; ibody++)
1687
 
    imagebody[ibody] = ((imageint) IMGMAX << IMG2BITS) | 
1688
 
      ((imageint) IMGMAX << IMGBITS) | IMGMAX;
 
1669
  if (!staticflag) {
 
1670
    for (ibody = 0; ibody < nbody; ibody++)
 
1671
      imagebody[ibody] = ((imageint) IMGMAX << IMG2BITS) | 
 
1672
        ((imageint) IMGMAX << IMGBITS) | IMGMAX;
 
1673
  }
1689
1674
 
1690
1675
  pre_neighbor();
1691
1676
 
1702
1687
    if (body[i] < 0) continue;
1703
1688
    ibody = body[i];
1704
1689
 
1705
 
    xbox = (image[i] & IMGMASK) - IMGMAX;
1706
 
    ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
1707
 
    zbox = (image[i] >> IMG2BITS) - IMGMAX;
 
1690
    xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
 
1691
    ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
 
1692
    zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
1708
1693
 
1709
1694
    if (triclinic == 0) {
1710
1695
      xunwrap = x[i][0] + xbox*xprd;
1857
1842
 
1858
1843
    ibody = body[i];
1859
1844
 
1860
 
    xbox = (image[i] & IMGMASK) - IMGMAX;
1861
 
    ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
1862
 
    zbox = (image[i] >> IMG2BITS) - IMGMAX;
 
1845
    xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
 
1846
    ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
 
1847
    zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
1863
1848
 
1864
1849
    if (triclinic == 0) {
1865
1850
      xunwrap = x[i][0] + xbox*xprd;
2018
2003
  }
2019
2004
 
2020
2005
  if (infile) memory->destroy(inbody);
 
2006
 
 
2007
  // static properties have now been initialized once
 
2008
  // used to prevent re-initialization which would re-read infile
 
2009
 
 
2010
  staticflag = 1;
2021
2011
}
2022
2012
 
2023
2013
/* ----------------------------------------------------------------------
2024
 
   one-time initialization of dynamic rigid body attributes
2025
 
   Vcm and angmom, computed explicitly from constituent particles
2026
 
   even if wrong for overlapping particles, is OK,
2027
 
     since is just setting initial time=0 Vcm and angmom of the body
 
2014
   initialization of dynamic rigid body attributes
 
2015
   set vcm and angmom, computed explicitly from constituent particles
 
2016
   OK if wrong for overlapping particles,
 
2017
     since is just setting vcm/angmom at start of run,
2028
2018
     which can be estimated value
2029
2019
------------------------------------------------------------------------- */
2030
2020
 
2040
2030
  double **v = atom->v;
2041
2031
  double *rmass = atom->rmass;
2042
2032
  double *mass = atom->mass;
2043
 
  imageint *image = atom->image;
2044
2033
  int *type = atom->type;
2045
2034
  int nlocal = atom->nlocal;
2046
2035
 
2061
2050
    sum[ibody][1] += v[i][1] * massone;
2062
2051
    sum[ibody][2] += v[i][2] * massone;
2063
2052
 
2064
 
    domain->unmap(x[i],image[i],unwrap);
 
2053
    domain->unmap(x[i],xcmimage[i],unwrap);
2065
2054
    dx = unwrap[0] - xcm[ibody][0];
2066
2055
    dy = unwrap[1] - xcm[ibody][1];
2067
2056
    dz = unwrap[2] - xcm[ibody][2];
2179
2168
    // id = rigid body ID
2180
2169
    // use ID as-is for SINGLE, as mol-ID for MOLECULE, as-is for GROUP
2181
2170
    // for which = 0, store mass/com in vec/array
2182
 
    // for which = 1, store interia tensor array, invert 3,4,5 values to Voigt
 
2171
    // for which = 1, store inertia tensor array, invert 3,4,5 values to Voigt
2183
2172
 
2184
2173
    for (int i = 0; i < nchunk; i++) {
2185
2174
      next = strchr(buf,'\n');
2282
2271
{
2283
2272
  int nmax = atom->nmax;
2284
2273
  double bytes = nmax * sizeof(int);
 
2274
  bytes += nmax * sizeof(imageint);
2285
2275
  bytes += nmax*3 * sizeof(double);
2286
2276
  bytes += maxvatom*6 * sizeof(double);    // vatom
2287
2277
  if (extended) {
2299
2289
void FixRigid::grow_arrays(int nmax)
2300
2290
{
2301
2291
  memory->grow(body,nmax,"rigid:body");
 
2292
  memory->grow(xcmimage,nmax,"rigid:xcmimage");
2302
2293
  memory->grow(displace,nmax,3,"rigid:displace");
2303
2294
  if (extended) {
2304
2295
    memory->grow(eflags,nmax,"rigid:eflags");
2314
2305
void FixRigid::copy_arrays(int i, int j, int delflag)
2315
2306
{
2316
2307
  body[j] = body[i];
 
2308
  xcmimage[j] = xcmimage[i];
2317
2309
  displace[j][0] = displace[i][0];
2318
2310
  displace[j][1] = displace[i][1];
2319
2311
  displace[j][2] = displace[i][2];
2336
2328
void FixRigid::set_arrays(int i)
2337
2329
{
2338
2330
  body[i] = -1;
 
2331
  xcmimage[i] = 0;
2339
2332
  displace[i][0] = 0.0;
2340
2333
  displace[i][1] = 0.0;
2341
2334
  displace[i][2] = 0.0;
2347
2340
 
2348
2341
int FixRigid::pack_exchange(int i, double *buf)
2349
2342
{
2350
 
  buf[0] = body[i];
2351
 
  buf[1] = displace[i][0];
2352
 
  buf[2] = displace[i][1];
2353
 
  buf[3] = displace[i][2];
2354
 
  if (!extended) return 4;
 
2343
  buf[0] = ubuf(body[i]).d;
 
2344
  buf[1] = ubuf(xcmimage[i]).d;
 
2345
  buf[2] = displace[i][0];
 
2346
  buf[3] = displace[i][1];
 
2347
  buf[4] = displace[i][2];
 
2348
  if (!extended) return 5;
2355
2349
 
2356
 
  int m = 4;
 
2350
  int m = 5;
2357
2351
  buf[m++] = eflags[i];
2358
2352
  for (int j = 0; j < orientflag; j++)
2359
2353
    buf[m++] = orient[i][j];
2371
2365
 
2372
2366
int FixRigid::unpack_exchange(int nlocal, double *buf)
2373
2367
{
2374
 
  body[nlocal] = static_cast<int> (buf[0]);
2375
 
  displace[nlocal][0] = buf[1];
2376
 
  displace[nlocal][1] = buf[2];
2377
 
  displace[nlocal][2] = buf[3];
2378
 
  if (!extended) return 4;
 
2368
  body[nlocal] = (int) ubuf(buf[0]).i;
 
2369
  xcmimage[nlocal] = (imageint) ubuf(buf[1]).i;
 
2370
  displace[nlocal][0] = buf[2];
 
2371
  displace[nlocal][1] = buf[3];
 
2372
  displace[nlocal][2] = buf[4];
 
2373
  if (!extended) return 5;
2379
2374
 
2380
 
  int m = 4;
 
2375
  int m = 5;
2381
2376
  eflags[nlocal] = static_cast<int> (buf[m++]);
2382
2377
  for (int j = 0; j < orientflag; j++)
2383
2378
    orient[nlocal][j] = buf[m++];