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