1 /* ----------------------------------------------------------------------
2     This is the
3 
4     ██╗     ██╗ ██████╗  ██████╗  ██████╗ ██╗  ██╗████████╗███████╗
5     ██║     ██║██╔════╝ ██╔════╝ ██╔════╝ ██║  ██║╚══██╔══╝██╔════╝
6     ██║     ██║██║  ███╗██║  ███╗██║  ███╗███████║   ██║   ███████╗
7     ██║     ██║██║   ██║██║   ██║██║   ██║██╔══██║   ██║   ╚════██║
8     ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║  ██║   ██║   ███████║
9     ╚══════╝╚═╝ ╚═════╝  ╚═════╝  ╚═════╝ ╚═╝  ╚═╝   ╚═╝   ╚══════╝®
10 
11     DEM simulation engine, released by
12     DCS Computing Gmbh, Linz, Austria
13     http://www.dcs-computing.com, office@dcs-computing.com
14 
15     LIGGGHTS® is part of CFDEM®project:
16     http://www.liggghts.com | http://www.cfdem.com
17 
18     Core developer and main author:
19     Christoph Kloss, christoph.kloss@dcs-computing.com
20 
21     LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22     License, version 2 or later. It is distributed in the hope that it will
23     be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25     received a copy of the GNU General Public License along with LIGGGHTS®.
26     If not, see http://www.gnu.org/licenses . See also top-level README
27     and LICENSE files.
28 
29     LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30     the producer of the LIGGGHTS® software and the CFDEM®coupling software
31     See http://www.cfdem.com/terms-trademark-policy for details.
32 
33 -------------------------------------------------------------------------
34     Contributing author and copyright for this file:
35     This file is from LAMMPS
36     LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37     http://lammps.sandia.gov, Sandia National Laboratories
38     Steve Plimpton, sjplimp@sandia.gov
39 
40     Copyright (2003) Sandia Corporation.  Under the terms of Contract
41     DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42     certain rights in this software.  This software is distributed under
43     the GNU General Public License.
44 ------------------------------------------------------------------------- */
45 
46 #include <stdlib.h>
47 #include "fix_minimize.h"
48 #include "atom.h"
49 #include "domain.h"
50 #include "memory.h"
51 #include "error.h"
52 
53 using namespace LAMMPS_NS;
54 using namespace FixConst;
55 
56 /* ---------------------------------------------------------------------- */
57 
FixMinimize(LAMMPS * lmp,int narg,char ** arg)58 FixMinimize::FixMinimize(LAMMPS *lmp, int narg, char **arg) :
59   Fix(lmp, narg, arg)
60 {
61   nvector = 0;
62   peratom = NULL;
63   vectors = NULL;
64 
65   // register callback to this fix from Atom class
66   // don't perform initial allocation here, must wait until add_vector()
67 
68   atom->add_callback(0);
69 }
70 
71 /* ---------------------------------------------------------------------- */
72 
~FixMinimize()73 FixMinimize::~FixMinimize()
74 {
75   // unregister callbacks to this fix from Atom class
76 
77   atom->delete_callback(id,0);
78 
79   // delete locally stored data
80 
81   memory->destroy(peratom);
82   for (int m = 0; m < nvector; m++) memory->destroy(vectors[m]);
83   memory->sfree(vectors);
84 }
85 
86 /* ---------------------------------------------------------------------- */
87 
setmask()88 int FixMinimize::setmask()
89 {
90   return 0;
91 }
92 
93 /* ----------------------------------------------------------------------
94    allocate/initialize memory for a new vector with N elements per atom
95 ------------------------------------------------------------------------- */
96 
add_vector(int n)97 void FixMinimize::add_vector(int n)
98 {
99   memory->grow(peratom,nvector+1,"minimize:peratom");
100   peratom[nvector] = n;
101 
102   vectors = (double **)
103     memory->srealloc(vectors,(nvector+1)*sizeof(double *),"minimize:vectors");
104   memory->create(vectors[nvector],atom->nmax*n,"minimize:vector");
105 
106   int ntotal = n*atom->nlocal;
107   for (int i = 0; i < ntotal; i++) vectors[nvector][i] = 0.0;
108   nvector++;
109 }
110 
111 /* ----------------------------------------------------------------------
112    return a pointer to the Mth vector
113 ------------------------------------------------------------------------- */
114 
request_vector(int m)115 double *FixMinimize::request_vector(int m)
116 {
117   return vectors[m];
118 }
119 
120 /* ----------------------------------------------------------------------
121    store box size at beginning of line search
122 ------------------------------------------------------------------------- */
123 
store_box()124 void FixMinimize::store_box()
125 {
126   boxlo[0] = domain->boxlo[0];
127   boxlo[1] = domain->boxlo[1];
128   boxlo[2] = domain->boxlo[2];
129   boxhi[0] = domain->boxhi[0];
130   boxhi[1] = domain->boxhi[1];
131   boxhi[2] = domain->boxhi[2];
132 }
133 
134 /* ----------------------------------------------------------------------
135    reset x0 for atoms that moved across PBC via reneighboring in line search
136    x0 = 1st vector
137    must do minimum_image using original box stored at beginning of line search
138    swap & set_global_box() change to original box, then restore current box
139 ------------------------------------------------------------------------- */
140 
reset_coords()141 void FixMinimize::reset_coords()
142 {
143   box_swap();
144   domain->set_global_box();
145 
146   double **x = atom->x;
147   double *x0 = vectors[0];
148   int nlocal = atom->nlocal;
149   double dx,dy,dz,dx0,dy0,dz0;
150 
151   int n = 0;
152   for (int i = 0; i < nlocal; i++) {
153     dx = dx0 = x[i][0] - x0[n];
154     dy = dy0 = x[i][1] - x0[n+1];
155     dz = dz0 = x[i][2] - x0[n+2];
156     domain->minimum_image(dx,dy,dz);
157     if (dx != dx0) x0[n] = x[i][0] - dx;
158     if (dy != dy0) x0[n+1] = x[i][1] - dy;
159     if (dz != dz0) x0[n+2] = x[i][2] - dz;
160     n += 3;
161   }
162 
163   box_swap();
164   domain->set_global_box();
165 }
166 
167 /* ----------------------------------------------------------------------
168    swap current box size with stored box size
169 ------------------------------------------------------------------------- */
170 
box_swap()171 void FixMinimize::box_swap()
172 {
173   double tmp;
174 
175   tmp = boxlo[0];
176   boxlo[0] = domain->boxlo[0];
177   domain->boxlo[0] = tmp;
178   tmp = boxlo[1];
179   boxlo[1] = domain->boxlo[1];
180   domain->boxlo[1] = tmp;
181   tmp = boxlo[2];
182   boxlo[2] = domain->boxlo[2];
183   domain->boxlo[2] = tmp;
184 
185   tmp = boxhi[0];
186   boxhi[0] = domain->boxhi[0];
187   domain->boxhi[0] = tmp;
188   tmp = boxhi[1];
189   boxhi[1] = domain->boxhi[1];
190   domain->boxhi[1] = tmp;
191   tmp = boxhi[2];
192   boxhi[2] = domain->boxhi[2];
193   domain->boxhi[2] = tmp;
194 }
195 
196 /* ----------------------------------------------------------------------
197    memory usage of local atom-based arrays
198 ------------------------------------------------------------------------- */
199 
memory_usage()200 double FixMinimize::memory_usage()
201 {
202   double bytes = 0.0;
203   for (int m = 0; m < nvector; m++)
204     bytes += atom->nmax*peratom[m]*sizeof(double);
205   return bytes;
206 }
207 
208 /* ----------------------------------------------------------------------
209    allocate local atom-based arrays
210 ------------------------------------------------------------------------- */
211 
grow_arrays(int nmax)212 void FixMinimize::grow_arrays(int nmax)
213 {
214   for (int m = 0; m < nvector; m++)
215     memory->grow(vectors[m],peratom[m]*nmax,"minimize:vector");
216 }
217 
218 /* ----------------------------------------------------------------------
219    copy values within local atom-based arrays
220 ------------------------------------------------------------------------- */
221 
copy_arrays(int i,int j,int delflag)222 void FixMinimize::copy_arrays(int i, int j, int delflag)
223 {
224   int m,iper,nper,ni,nj;
225 
226   for (m = 0; m < nvector; m++) {
227     nper = peratom[m];
228     ni = nper*i;
229     nj = nper*j;
230     for (iper = 0; iper < nper; iper++) vectors[m][nj++] = vectors[m][ni++];
231   }
232 }
233 
234 /* ----------------------------------------------------------------------
235    pack values in local atom-based arrays for exchange with another proc
236 ------------------------------------------------------------------------- */
237 
pack_exchange(int i,double * buf)238 int FixMinimize::pack_exchange(int i, double *buf)
239 {
240   int m,iper,nper,ni;
241 
242   int n = 0;
243   for (m = 0; m < nvector; m++) {
244     nper = peratom[m];
245     ni = nper*i;
246     for (iper = 0; iper < nper; iper++) buf[n++] = vectors[m][ni++];
247   }
248   return n;
249 }
250 
251 /* ----------------------------------------------------------------------
252    unpack values in local atom-based arrays from exchange with another proc
253 ------------------------------------------------------------------------- */
254 
unpack_exchange(int nlocal,double * buf)255 int FixMinimize::unpack_exchange(int nlocal, double *buf)
256 {
257   int m,iper,nper,ni;
258 
259   int n = 0;
260   for (m = 0; m < nvector; m++) {
261     nper = peratom[m];
262     ni = nper*nlocal;
263     for (iper = 0; iper < nper; iper++) vectors[m][ni++] = buf[n++];
264   }
265   return n;
266 }
267