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 "delete_bonds.h"
16
17 #include "atom.h"
18 #include "atom_vec.h"
19 #include "comm.h"
20 #include "domain.h"
21 #include "error.h"
22 #include "force.h"
23 #include "group.h"
24 #include "special.h"
25
26 #include <cstring>
27
28 using namespace LAMMPS_NS;
29
30 enum{MULTI,ATOM,BOND,ANGLE,DIHEDRAL,IMPROPER,STATS};
31
32 /* ---------------------------------------------------------------------- */
33
DeleteBonds(LAMMPS * lmp)34 DeleteBonds::DeleteBonds(LAMMPS *lmp) : Command(lmp) {}
35
36 /* ---------------------------------------------------------------------- */
37
command(int narg,char ** arg)38 void DeleteBonds::command(int narg, char **arg)
39 {
40 if (domain->box_exist == 0)
41 error->all(FLERR,"Delete_bonds command before simulation box is defined");
42 if (atom->natoms == 0)
43 error->all(FLERR,"Delete_bonds command with no atoms existing");
44 if (atom->molecular != Atom::MOLECULAR)
45 error->all(FLERR,"Cannot use delete_bonds with non-molecular system");
46
47 if (narg < 2) error->all(FLERR,"Illegal delete_bonds command");
48
49 // init entire system since comm->borders is done
50 // comm::init needs neighbor::init needs pair::init needs kspace::init, etc
51
52 if (comm->me == 0) utils::logmesg(lmp,"System init for delete_bonds ...\n");
53 lmp->init();
54
55 if (comm->me == 0) utils::logmesg(lmp,"Deleting bonds ...\n");
56
57 // identify group
58
59 int igroup = group->find(arg[0]);
60 if (igroup == -1) error->all(FLERR,"Cannot find delete_bonds group ID");
61 int groupbit = group->bitmask[igroup];
62
63 // set style and which = type value
64
65 int style = -1;
66 if (strcmp(arg[1],"multi") == 0) style = MULTI;
67 else if (strcmp(arg[1],"atom") == 0) style = ATOM;
68 else if (strcmp(arg[1],"bond") == 0) style = BOND;
69 else if (strcmp(arg[1],"angle") == 0) style = ANGLE;
70 else if (strcmp(arg[1],"dihedral") == 0) style = DIHEDRAL;
71 else if (strcmp(arg[1],"improper") == 0) style = IMPROPER;
72 else if (strcmp(arg[1],"stats") == 0) style = STATS;
73 else error->all(FLERR,"Illegal delete_bonds command");
74
75 // setup list of types (atom,bond,etc) to consider
76 // use utils::bounds(FLERR,) to allow setting of range of types
77 // range can be 0 to ntypes inclusive
78
79 int *tlist = nullptr;
80
81 int iarg = 2;
82 if (style != MULTI && style != STATS) {
83 if (narg < 3) error->all(FLERR,"Illegal delete_bonds command");
84
85 int n = -1;
86 if (style == ATOM) n = atom->ntypes;
87 if (style == BOND) n = atom->nbondtypes;
88 if (style == ANGLE) n = atom->nangletypes;
89 if (style == DIHEDRAL) n = atom->ndihedraltypes;
90 if (style == IMPROPER) n = atom->nimpropertypes;
91
92 tlist = new int[n+1];
93 for (int i = 0; i <= n; i++) tlist[i] = 0;
94 int nlo,nhi;
95 utils::bounds(FLERR,arg[2],0,n,nlo,nhi,error);
96 for (int i = nlo; i <= nhi; i++) tlist[i] = 1;
97
98 iarg++;
99 }
100
101 // grab optional keywords
102
103 int any_flag = 0;
104 int undo_flag = 0;
105 int remove_flag = 0;
106 int special_flag = 0;
107 int induce_flag = 0;
108
109 while (iarg < narg) {
110 if (strcmp(arg[iarg],"any") == 0) any_flag = 1;
111 else if (strcmp(arg[iarg],"undo") == 0) undo_flag = 1;
112 else if (strcmp(arg[iarg],"remove") == 0) remove_flag = 1;
113 else if (strcmp(arg[iarg],"special") == 0) special_flag = 1;
114 else if (strcmp(arg[iarg],"induce") == 0) induce_flag = 1;
115 else error->all(FLERR,"Illegal delete_bonds command");
116 iarg++;
117 }
118
119 // border swap to insure type and mask is current for off-proc atoms
120 // enforce PBC before in case atoms are outside box
121
122 if (domain->triclinic) domain->x2lamda(atom->nlocal);
123 domain->pbc();
124 domain->reset_box();
125 comm->setup();
126 comm->exchange();
127 comm->borders();
128 if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
129
130 // set topology interactions either off or on
131 // criteria for an interaction to potentially be changed (set flag = 1)
132 // all atoms or any atom in interaction must be in group, based on any_flag
133 // for style = MULTI, all bond/angle/dihedral/improper, no other criteria
134 // for style = ATOM, same as MULTI, plus at least one atom is specified type
135 // for style = BOND/ANGLE/DIHEDRAL/IMPROPER, interaction is specified type
136 // for style = STATS only compute stats, flag is always 0
137 // if flag = 1
138 // set interaction type negative if undo_flag = 0
139 // set interaction type positive if undo_flag = 1
140
141 int *mask = atom->mask;
142 int *type = atom->type;
143 int nlocal = atom->nlocal;
144
145 int i,m,n,consider,flag,itype;
146 int atom1,atom2,atom3,atom4;
147
148 if (atom->avec->bonds_allow &&
149 (style == BOND || style == MULTI || style == ATOM)) {
150 int *num_bond = atom->num_bond;
151 int **bond_type = atom->bond_type;
152
153 for (i = 0; i < nlocal; i++) {
154 for (m = 0; m < num_bond[i]; m++) {
155 atom1 = atom->map(atom->bond_atom[i][m]);
156 if (atom1 == -1) error->one(FLERR,"Bond atom missing in delete_bonds");
157 consider = 0;
158 if (!any_flag && mask[i] & groupbit && mask[atom1] & groupbit)
159 consider = 1;
160 if (any_flag && (mask[i] & groupbit || mask[atom1] & groupbit))
161 consider = 1;
162 if (consider) {
163 flag = 0;
164 if (style == MULTI) flag = 1;
165 else if (style == ATOM) {
166 if (tlist[type[i]] || tlist[type[atom1]]) flag = 1;
167 } else if (style == BOND) {
168 itype = abs(bond_type[i][m]);
169 if (tlist[itype]) flag = 1;
170 }
171 if (flag) {
172 if (undo_flag == 0 && bond_type[i][m] > 0)
173 bond_type[i][m] = -bond_type[i][m];
174 if (undo_flag == 1 && bond_type[i][m] < 0)
175 bond_type[i][m] = -bond_type[i][m];
176 }
177 }
178 }
179 }
180 }
181
182 if (atom->avec->angles_allow &&
183 (style == ANGLE || style == MULTI || style == ATOM)) {
184 int *num_angle = atom->num_angle;
185 int **angle_type = atom->angle_type;
186
187 for (i = 0; i < nlocal; i++) {
188 for (m = 0; m < num_angle[i]; m++) {
189 atom1 = atom->map(atom->angle_atom1[i][m]);
190 atom2 = atom->map(atom->angle_atom2[i][m]);
191 atom3 = atom->map(atom->angle_atom3[i][m]);
192 if (atom1 == -1 || atom2 == -1 || atom3 == -1)
193 error->one(FLERR,"Angle atom missing in delete_bonds");
194 consider = 0;
195 if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
196 mask[atom3] & groupbit) consider = 1;
197 if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
198 mask[atom3] & groupbit)) consider = 1;
199 if (consider) {
200 flag = 0;
201 if (style == MULTI) flag = 1;
202 else if (style == ATOM) {
203 if (tlist[type[atom1]] || tlist[type[atom2]] ||
204 tlist[type[atom3]]) flag = 1;
205 } else if (style == ANGLE) {
206 itype = abs(angle_type[i][m]);
207 if (tlist[itype]) flag = 1;
208 }
209 if (flag) {
210 if (undo_flag == 0 && angle_type[i][m] > 0)
211 angle_type[i][m] = -angle_type[i][m];
212 if (undo_flag == 1 && angle_type[i][m] < 0)
213 angle_type[i][m] = -angle_type[i][m];
214 }
215 }
216 }
217 }
218 }
219
220 if (atom->avec->dihedrals_allow &&
221 (style == DIHEDRAL || style == MULTI || style == ATOM)) {
222 int *num_dihedral = atom->num_dihedral;
223 int **dihedral_type = atom->dihedral_type;
224
225 for (i = 0; i < nlocal; i++) {
226 for (m = 0; m < num_dihedral[i]; m++) {
227 atom1 = atom->map(atom->dihedral_atom1[i][m]);
228 atom2 = atom->map(atom->dihedral_atom2[i][m]);
229 atom3 = atom->map(atom->dihedral_atom3[i][m]);
230 atom4 = atom->map(atom->dihedral_atom4[i][m]);
231 if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
232 error->one(FLERR,"Dihedral atom missing in delete_bonds");
233 consider = 0;
234 if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
235 mask[atom3] & groupbit && mask[atom4] & groupbit) consider = 1;
236 if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
237 mask[atom3] & groupbit || mask[atom4] & groupbit))
238 consider = 1;
239 if (consider) {
240 flag = 0;
241 if (style == MULTI) flag = 1;
242 else if (style == ATOM) {
243 if (tlist[type[atom1]] || tlist[type[atom2]] ||
244 tlist[type[atom3]] || tlist[type[atom4]]) flag = 1;
245 } else if (style == DIHEDRAL) {
246 itype = abs(dihedral_type[i][m]);
247 if (tlist[itype]) flag = 1;
248 }
249 if (flag) {
250 if (undo_flag == 0 && dihedral_type[i][m] > 0)
251 dihedral_type[i][m] = -dihedral_type[i][m];
252 if (undo_flag == 1 && dihedral_type[i][m] < 0)
253 dihedral_type[i][m] = -dihedral_type[i][m];
254 }
255 }
256 }
257 }
258 }
259
260 if (atom->avec->impropers_allow &&
261 (style == IMPROPER || style == MULTI || style == ATOM)) {
262 int *num_improper = atom->num_improper;
263 int **improper_type = atom->improper_type;
264
265 for (i = 0; i < nlocal; i++) {
266 for (m = 0; m < num_improper[i]; m++) {
267 atom1 = atom->map(atom->improper_atom1[i][m]);
268 atom2 = atom->map(atom->improper_atom2[i][m]);
269 atom3 = atom->map(atom->improper_atom3[i][m]);
270 atom4 = atom->map(atom->improper_atom4[i][m]);
271 if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
272 error->one(FLERR,"Improper atom missing in delete_bonds");
273 consider = 0;
274 if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
275 mask[atom3] & groupbit && mask[atom4] & groupbit) consider = 1;
276 if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
277 mask[atom3] & groupbit || mask[atom4] & groupbit))
278 consider = 1;
279 if (consider) {
280 flag = 0;
281 if (style == MULTI) flag = 1;
282 else if (style == ATOM) {
283 if (tlist[type[atom1]] || tlist[type[atom2]] ||
284 tlist[type[atom3]] || tlist[type[atom4]]) flag = 1;
285 } else if (style == IMPROPER) {
286 itype = abs(improper_type[i][m]);
287 if (tlist[itype]) flag = 1;
288 }
289 if (flag) {
290 if (undo_flag == 0 && improper_type[i][m] > 0)
291 improper_type[i][m] = -improper_type[i][m];
292 if (undo_flag == 1 && improper_type[i][m] < 0)
293 improper_type[i][m] = -improper_type[i][m];
294 }
295 }
296 }
297 }
298 }
299
300 delete [] tlist;
301
302 // induce turn off of angles, dihedral, impropers due to turned off bonds
303 // induce turn off of dihedrals due to turned off angles
304 // all atoms or any atom in interaction must be in group, based on any_flag
305
306 if (induce_flag) {
307
308 // circulate list of turned off bonds around ring of procs
309
310 // circulate list of turned off angles around ring of procs
311
312 }
313
314 // remove interactions if requested
315 // all atoms or any atom in interaction must be in group, based on any_flag
316
317 if (remove_flag) {
318
319 if (atom->avec->bonds_allow) {
320 for (i = 0; i < nlocal; i++) {
321 m = 0;
322 while (m < atom->num_bond[i]) {
323 if (atom->bond_type[i][m] <= 0) {
324 atom1 = atom->map(atom->bond_atom[i][m]);
325 flag = 0;
326 if (!any_flag && mask[i] & groupbit && mask[atom1] & groupbit)
327 flag = 1;
328 if (any_flag && (mask[i] & groupbit || mask[atom1] & groupbit))
329 flag = 1;
330 if (flag) {
331 n = atom->num_bond[i];
332 atom->bond_type[i][m] = atom->bond_type[i][n-1];
333 atom->bond_atom[i][m] = atom->bond_atom[i][n-1];
334 atom->num_bond[i]--;
335 } else m++;
336 } else m++;
337 }
338 }
339 }
340
341 if (atom->avec->angles_allow) {
342 for (i = 0; i < nlocal; i++) {
343 m = 0;
344 while (m < atom->num_angle[i]) {
345 if (atom->angle_type[i][m] <= 0) {
346 atom1 = atom->map(atom->angle_atom1[i][m]);
347 atom2 = atom->map(atom->angle_atom2[i][m]);
348 atom3 = atom->map(atom->angle_atom3[i][m]);
349 flag = 0;
350 if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
351 mask[atom3] & groupbit) flag = 1;
352 if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
353 mask[atom3] & groupbit)) flag = 1;
354 if (flag) {
355 n = atom->num_angle[i];
356 atom->angle_type[i][m] = atom->angle_type[i][n-1];
357 atom->angle_atom1[i][m] = atom->angle_atom1[i][n-1];
358 atom->angle_atom2[i][m] = atom->angle_atom2[i][n-1];
359 atom->angle_atom3[i][m] = atom->angle_atom3[i][n-1];
360 atom->num_angle[i]--;
361 } else m++;
362 } else m++;
363 }
364 }
365 }
366
367 if (atom->avec->dihedrals_allow) {
368 for (i = 0; i < nlocal; i++) {
369 m = 0;
370 while (m < atom->num_dihedral[i]) {
371 if (atom->dihedral_type[i][m] <= 0) {
372 atom1 = atom->map(atom->dihedral_atom1[i][m]);
373 atom2 = atom->map(atom->dihedral_atom2[i][m]);
374 atom3 = atom->map(atom->dihedral_atom3[i][m]);
375 atom4 = atom->map(atom->dihedral_atom4[i][m]);
376 flag = 0;
377 if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
378 mask[atom3] & groupbit && mask[atom4] & groupbit) flag = 1;
379 if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
380 mask[atom3] & groupbit || mask[atom4] & groupbit))
381 flag = 1;
382 if (flag) {
383 n = atom->num_dihedral[i];
384 atom->dihedral_type[i][m] = atom->dihedral_type[i][n-1];
385 atom->dihedral_atom1[i][m] = atom->dihedral_atom1[i][n-1];
386 atom->dihedral_atom2[i][m] = atom->dihedral_atom2[i][n-1];
387 atom->dihedral_atom3[i][m] = atom->dihedral_atom3[i][n-1];
388 atom->dihedral_atom4[i][m] = atom->dihedral_atom4[i][n-1];
389 atom->num_dihedral[i]--;
390 } else m++;
391 } else m++;
392 }
393 }
394 }
395
396 if (atom->avec->impropers_allow) {
397 for (i = 0; i < nlocal; i++) {
398 m = 0;
399 while (m < atom->num_improper[i]) {
400 if (atom->improper_type[i][m] <= 0) {
401 atom1 = atom->map(atom->improper_atom1[i][m]);
402 atom2 = atom->map(atom->improper_atom2[i][m]);
403 atom3 = atom->map(atom->improper_atom3[i][m]);
404 atom4 = atom->map(atom->improper_atom4[i][m]);
405 flag = 0;
406 if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
407 mask[atom3] & groupbit && mask[atom4] & groupbit) flag = 1;
408 if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
409 mask[atom3] & groupbit || mask[atom4] & groupbit))
410 flag = 1;
411 if (flag) {
412 n = atom->num_improper[i];
413 atom->improper_type[i][m] = atom->improper_type[i][n-1];
414 atom->improper_atom1[i][m] = atom->improper_atom1[i][n-1];
415 atom->improper_atom2[i][m] = atom->improper_atom2[i][n-1];
416 atom->improper_atom3[i][m] = atom->improper_atom3[i][n-1];
417 atom->improper_atom4[i][m] = atom->improper_atom4[i][n-1];
418 atom->num_improper[i]--;
419 } else m++;
420 } else m++;
421 }
422 }
423 }
424
425 }
426
427 // if interactions were removed, recompute global counts
428
429 if (remove_flag) {
430
431 if (atom->avec->bonds_allow) {
432 bigint nbonds = 0;
433 for (i = 0; i < nlocal; i++) nbonds += atom->num_bond[i];
434 MPI_Allreduce(&nbonds,&atom->nbonds,1,MPI_LMP_BIGINT,
435 MPI_SUM,world);
436 if (force->newton_bond == 0) atom->nbonds /= 2;
437 }
438
439 if (atom->avec->angles_allow) {
440 bigint nangles = 0;
441 for (i = 0; i < nlocal; i++) nangles += atom->num_angle[i];
442 MPI_Allreduce(&nangles,&atom->nangles,1,MPI_LMP_BIGINT,
443 MPI_SUM,world);
444 if (force->newton_bond == 0) atom->nangles /= 3;
445 }
446
447 if (atom->avec->dihedrals_allow) {
448 bigint ndihedrals = 0;
449 for (i = 0; i < nlocal; i++) ndihedrals += atom->num_dihedral[i];
450 MPI_Allreduce(&ndihedrals,&atom->ndihedrals,
451 1,MPI_LMP_BIGINT,MPI_SUM,world);
452 if (force->newton_bond == 0) atom->ndihedrals /= 4;
453 }
454
455 if (atom->avec->impropers_allow) {
456 bigint nimpropers = 0;
457 for (i = 0; i < nlocal; i++) nimpropers += atom->num_improper[i];
458 MPI_Allreduce(&nimpropers,&atom->nimpropers,
459 1,MPI_LMP_BIGINT,MPI_SUM,world);
460 if (force->newton_bond == 0) atom->nimpropers /= 4;
461 }
462
463 }
464
465 // compute and print stats
466
467 bigint tmp;
468 bigint bond_on,bond_off;
469 bigint angle_on,angle_off;
470 bigint dihedral_on,dihedral_off;
471 bigint improper_on,improper_off;
472
473 if (atom->avec->bonds_allow) {
474 bond_on = bond_off = 0;
475 for (i = 0; i < nlocal; i++)
476 for (m = 0; m < atom->num_bond[i]; m++)
477 if (atom->bond_type[i][m] > 0) bond_on++;
478 else bond_off++;
479 MPI_Allreduce(&bond_on,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
480 bond_on = tmp;
481 MPI_Allreduce(&bond_off,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
482 bond_off = tmp;
483 if (force->newton_bond == 0) {
484 bond_on /= 2;
485 bond_off /= 2;
486 }
487 }
488
489 if (atom->avec->angles_allow) {
490 angle_on = angle_off = 0;
491 for (i = 0; i < nlocal; i++)
492 for (m = 0; m < atom->num_angle[i]; m++)
493 if (atom->angle_type[i][m] > 0) angle_on++;
494 else angle_off++;
495 MPI_Allreduce(&angle_on,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
496 angle_on = tmp;
497 MPI_Allreduce(&angle_off,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
498 angle_off = tmp;
499 if (force->newton_bond == 0) {
500 angle_on /= 3;
501 angle_off /= 3;
502 }
503 }
504
505 if (atom->avec->dihedrals_allow) {
506 dihedral_on = dihedral_off = 0;
507 for (i = 0; i < nlocal; i++)
508 for (m = 0; m < atom->num_dihedral[i]; m++)
509 if (atom->dihedral_type[i][m] > 0) dihedral_on++;
510 else dihedral_off++;
511 MPI_Allreduce(&dihedral_on,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
512 dihedral_on = tmp;
513 MPI_Allreduce(&dihedral_off,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
514 dihedral_off = tmp;
515 if (force->newton_bond == 0) {
516 dihedral_on /= 4;
517 dihedral_off /= 4;
518 }
519 }
520
521 if (atom->avec->impropers_allow) {
522 improper_on = improper_off = 0;
523 for (i = 0; i < nlocal; i++)
524 for (m = 0; m < atom->num_improper[i]; m++)
525 if (atom->improper_type[i][m] > 0) improper_on++;
526 else improper_off++;
527 MPI_Allreduce(&improper_on,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
528 improper_on = tmp;
529 MPI_Allreduce(&improper_off,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
530 improper_off = tmp;
531 if (force->newton_bond == 0) {
532 improper_on /= 4;
533 improper_off /= 4;
534 }
535 }
536
537 if (comm->me == 0) {
538 if (atom->avec->bonds_allow)
539 utils::logmesg(lmp," {} total bonds, {} turned on, {} turned off\n",
540 atom->nbonds,bond_on,bond_off);
541
542 if (atom->avec->angles_allow)
543 utils::logmesg(lmp," {} total angles, {} turned on, {} turned off\n",
544 atom->nangles,angle_on,angle_off);
545
546 if (atom->avec->dihedrals_allow)
547 utils::logmesg(lmp," {} total dihedrals, {} turned on, {} turned off\n",
548 atom->ndihedrals,dihedral_on,dihedral_off);
549
550 if (atom->avec->impropers_allow)
551 utils::logmesg(lmp," {} total impropers, {} turned on, {} turned off\n",
552 atom->nimpropers,improper_on,improper_off);
553 }
554
555 // re-compute special list if requested
556
557 if (special_flag) {
558 Special special(lmp);
559 special.build();
560 }
561 }
562