← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/atc/ATC_HardyKernel.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2013-11-20 22:41:36 UTC
  • mfrom: (1.2.2)
  • Revision ID: package-import@ubuntu.com-20131120224136-tzx7leh606fqnckm
Tags: 0~20131119.git7162cf0-1
* [e65b919] Imported Upstream version 0~20131119.git7162cf0
* [f7bddd4] Fix some problems, introduced by upstream recently.
* [3616dfc] Use wrap-and-sort script.
* [7e92030] Ignore quilt dir

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/atc/ATC_HardyKernel.cpp
1
 
#include "ATC_HardyKernel.h"
2
 
#include "math.h"
3
 
#include <iostream> //for debugging purposes; take this out after I'm done
4
 
#include <vector>
5
 
#include "ATC_Error.h"
6
 
#include "Quadrature.h"
7
 
 
8
 
using namespace std;
9
 
 
10
 
static const double Pi = 4.0*atan(1.0);
11
 
static const double tol = 1.0e-8;
12
 
 
13
 
namespace ATC {
14
 
 
15
 
 
16
 
  //------------------------------------------------------------------------
17
 
  // constructor
18
 
  ATC_HardyKernel::ATC_HardyKernel(int nparameters, double* parameters):
19
 
    lammpsInterface_(LammpsInterface::instance()),
20
 
    Rc_(0),invRc_(0),nsd_(3)
21
 
  { 
22
 
    Rc_ = parameters[0];
23
 
    invRc_ = 1.0/Rc_;
24
 
    Rc_ = parameters[0];
25
 
    invRc_ = 1.0/Rc_;
26
 
    invVol_ = 1.0/(4.0/3.0*Pi*pow(Rc_,3));
27
 
 
28
 
    set_line_quadrature(line_ngauss,line_xg,line_wg);
29
 
 
30
 
    // get periodicity and box bounds/lengths
31
 
    lammpsInterface_->get_box_periodicity(periodicity[0],
32
 
                                          periodicity[1],periodicity[2]);
33
 
    lammpsInterface_->get_box_bounds(box_bounds[0][0],box_bounds[1][0],
34
 
                                     box_bounds[0][1],box_bounds[1][1],
35
 
                                     box_bounds[0][2],box_bounds[1][2]);
36
 
    for (int k = 0; k < 3; k++) {
37
 
      box_length[k] = box_bounds[1][k] - box_bounds[0][k]; 
38
 
    } 
39
 
  };
40
 
 
41
 
  // bond function value via quadrature
42
 
  double ATC_HardyKernel::bond(DENS_VEC& xa, DENS_VEC&xb, double lam1, double lam2)
43
 
  {
44
 
    DENS_VEC xab(nsd_), q(nsd_);
45
 
    double lamg;
46
 
    double bhsum=0.0;
47
 
    xab = xa - xb;
48
 
    for (int i = 0; i < line_ngauss; i++) {
49
 
      lamg=0.5*((lam2-lam1)*line_xg[i]+(lam2+lam1));
50
 
      q = lamg*xab + xb;
51
 
      double locg_value=this->value(q);
52
 
      bhsum+=locg_value*line_wg[i];
53
 
    }
54
 
    return 0.5*(lam2-lam1)*bhsum;
55
 
  }
56
 
 
57
 
  // localization-volume intercepts for bond calculation 
58
 
  // bond intercept values assuming spherical support
59
 
  void ATC_HardyKernel::bond_intercepts(DENS_VEC& xa,
60
 
                 DENS_VEC& xb, double &lam1, double &lam2) 
61
 
  {
62
 
    if (nsd_ == 2) {// for cylinders, axis is always z! 
63
 
      const int iaxis = 2;
64
 
      xa[iaxis] = 0.0;
65
 
      xb[iaxis] = 0.0;
66
 
    }
67
 
    lam1 = lam2 = -1;
68
 
    double ra_n = invRc_*xa.norm(); // lambda = 1
69
 
    double rb_n = invRc_*xb.norm(); // lambda = 0
70
 
    bool a_in = (ra_n <= 1.0);
71
 
    bool b_in = (rb_n <= 1.0);
72
 
    if (a_in && b_in) {
73
 
      lam1 = 0.0; 
74
 
      lam2 = 1.0; 
75
 
      return;
76
 
    }
77
 
    DENS_VEC xab = xa - xb;
78
 
    double rab_n = invRc_*xab.norm();
79
 
    double a = rab_n*rab_n; // always at least an interatomic distance
80
 
    double b = 2.0*invRc_*invRc_*xab.dot(xb);
81
 
    double c = rb_n*rb_n - 1.0;
82
 
    double discrim = b*b - 4.0*a*c; // discriminant
83
 
    if (discrim < 0) return; // no intersection
84
 
    // num recipes:
85
 
    double s1, s2;
86
 
    if (b < 0.0) {
87
 
      double aux = -0.5*(b-sqrt(discrim));
88
 
      s1 = c/aux;
89
 
      s2 = aux/a;
90
 
    } 
91
 
    else {
92
 
      double aux = -0.5*(b+sqrt(discrim));
93
 
      s1 = aux/a;
94
 
      s2 = c/aux;
95
 
    }
96
 
    if (a_in && !b_in) {
97
 
      lam1 = s1;
98
 
      lam2 = 1.0;
99
 
    } 
100
 
    else if (!a_in && b_in) { 
101
 
      lam1 = 0.0;
102
 
      lam2 = s2; 
103
 
    }
104
 
    else {
105
 
      if (s1 >= 0.0 && s2 <= 1.0) { 
106
 
        lam1 = s1; 
107
 
        lam2 = s2; 
108
 
      }
109
 
    }
110
 
  };
111
 
 
112
 
  //------------------------------------------------------------------------
113
 
  // constructor
114
 
  ATC_HardyKernelStep::ATC_HardyKernelStep
115
 
    (int nparameters, double* parameters): 
116
 
    ATC_HardyKernel(nparameters, parameters) 
117
 
  {
118
 
    for (int k = 0; k < nsd_; k++ ) {
119
 
      if ((bool) periodicity[k]) {
120
 
        if (Rc_ > 0.5*box_length[k]) {
121
 
          throw ATC_Error(0,"Size of localization volume is too large for periodic boundary condition");
122
 
        };
123
 
      };
124
 
    }; 
125
 
  }
126
 
 
127
 
  // function value
128
 
  double ATC_HardyKernelStep::value(DENS_VEC& x_atom) 
129
 
  {
130
 
    double rn=invRc_*x_atom.norm();
131
 
    if (rn <= 1.0) { return 1.0; }
132
 
    else           { return 0.0; }
133
 
  };
134
 
 
135
 
  //------------------------------------------------------------------------
136
 
  /** a step with rectangular support suitable for a rectangular grid */
137
 
  // constructor
138
 
  ATC_HardyKernelCell::ATC_HardyKernelCell
139
 
    (int nparameters, double* parameters): 
140
 
    ATC_HardyKernel(nparameters, parameters) 
141
 
  {
142
 
    hx = parameters[0];
143
 
    hy = parameters[1];
144
 
    hz = parameters[2];
145
 
    invVol_ = 1.0/8.0/(hx*hy*hz);
146
 
    cellBounds_.reset(6);
147
 
    cellBounds_(0) = -hx;
148
 
    cellBounds_(1) =  hx;
149
 
    cellBounds_(2) = -hy;
150
 
    cellBounds_(3) =  hy;
151
 
    cellBounds_(4) = -hz;
152
 
    cellBounds_(5) =  hz;
153
 
      
154
 
    for (int k = 0; k < nsd_; k++ ) {
155
 
      if ((bool) periodicity[k]) {
156
 
        if (parameters[k] > 0.5*box_length[k]) {
157
 
          throw ATC_Error(0,"Size of localization volume is too large for periodic boundary condition");
158
 
        };
159
 
      };
160
 
    };
161
 
  }
162
 
 
163
 
  // function value
164
 
  double ATC_HardyKernelCell::value(DENS_VEC& x_atom) 
165
 
  {
166
 
    if ((cellBounds_(0) <= x_atom(0)) && (x_atom(0) < cellBounds_(1)) 
167
 
     && (cellBounds_(2) <= x_atom(1)) && (x_atom(1) < cellBounds_(3)) 
168
 
     && (cellBounds_(4) <= x_atom(2)) && (x_atom(2) < cellBounds_(5))) { 
169
 
      return 1.0; 
170
 
    } 
171
 
    else { 
172
 
      return 0.0; 
173
 
    }
174
 
  };
175
 
 
176
 
  // bond intercept values for rectangular region : origin is the node position
177
 
  void ATC_HardyKernelCell::bond_intercepts(DENS_VEC& xa,
178
 
                 DENS_VEC& xb, double &lam1, double &lam2) 
179
 
  {
180
 
    lam1 = 0.0; // start
181
 
    lam2 = 1.0; // end
182
 
 
183
 
    bool a_in = (value(xa) > 0.0);
184
 
    bool b_in = (value(xb) > 0.0);
185
 
 
186
 
    // (1) both in, no intersection needed
187
 
    if (a_in && b_in) {
188
 
      return;
189
 
    }
190
 
    // (2) for one in & one out -> one plane intersection
191
 
    // determine the points of intersection between the line joining
192
 
    // atoms a and b and the bounding planes of the localization volume
193
 
    else if (a_in || b_in) {
194
 
      DENS_VEC xab = xa - xb;
195
 
      for (int i = 0; i < nsd_; i++) {
196
 
        // check if segment is parallel to face
197
 
        if (fabs(xab(i)) > tol) {
198
 
          for (int j = 0; j < 2; j++) {
199
 
            double s = (cellBounds_(2*i+j) - xb(i))/xab(i);
200
 
            // check if between a & b
201
 
            if (s >= 0 && s <= 1) {
202
 
              bool in_bounds = false; 
203
 
              DENS_VEC x = xb + s*xab;
204
 
              if (i == 0) { 
205
 
                if ((cellBounds_(2) <= x(1)) && (x(1) <= cellBounds_(3))
206
 
                 && (cellBounds_(4) <= x(2)) && (x(2) <= cellBounds_(5))) {
207
 
                  in_bounds = true;
208
 
                }
209
 
              }
210
 
              else if (i == 1) { 
211
 
                if ((cellBounds_(0) <= x(0)) && (x(0) <= cellBounds_(1))
212
 
                 && (cellBounds_(4) <= x(2)) && (x(2) <= cellBounds_(5))) {
213
 
                  in_bounds = true;
214
 
                }
215
 
              }
216
 
              else if (i == 2) { 
217
 
                if ((cellBounds_(0) <= x(0)) && (x(0) <= cellBounds_(1))
218
 
                 && (cellBounds_(2) <= x(1)) && (x(1) <= cellBounds_(3))) {
219
 
                  in_bounds = true;
220
 
                }
221
 
              }
222
 
              if (in_bounds) {
223
 
                if (a_in) { lam1 = s;}
224
 
                else      { lam2 = s;}
225
 
                return;
226
 
              }
227
 
            }
228
 
          }
229
 
        } 
230
 
      }
231
 
      throw ATC_Error(0,"logic failure in HardyKernel Cell for single intersection\n");
232
 
    }
233
 
    // (3) both out -> corner intersection
234
 
    else {
235
 
      lam2 = lam1; // default to no intersection
236
 
      DENS_VEC xab = xa - xb;
237
 
      double ss[6] = {-1,-1,-1,-1,-1,-1};
238
 
      int is = 0;
239
 
      for (int i = 0; i < nsd_; i++) {
240
 
        // check if segment is parallel to face
241
 
        if (fabs(xab(i)) > tol) {
242
 
          for (int j = 0; j < 2; j++) {
243
 
            double s = (cellBounds_(2*i+j) - xb(i))/xab(i);
244
 
            // check if between a & b
245
 
            if (s >= 0 && s <= 1) {
246
 
              // check if in face
247
 
              DENS_VEC x = xb + s*xab;
248
 
              if (i == 0) { 
249
 
                if ((cellBounds_(2) <= x(1)) && (x(1) <= cellBounds_(3))
250
 
                 && (cellBounds_(4) <= x(2)) && (x(2) <= cellBounds_(5))) {
251
 
                  ss[is++] = s;
252
 
                }
253
 
              }
254
 
              else if (i == 1) { 
255
 
                if ((cellBounds_(0) <= x(0)) && (x(0) <= cellBounds_(1))
256
 
                 && (cellBounds_(4) <= x(2)) && (x(2) <= cellBounds_(5))) {
257
 
                  ss[is++] = s;
258
 
                }
259
 
              }
260
 
              else if (i == 2) { 
261
 
                if ((cellBounds_(0) <= x(0)) && (x(0) <= cellBounds_(1))
262
 
                 && (cellBounds_(2) <= x(1)) && (x(1) <= cellBounds_(3))) {
263
 
                  ss[is++] = s;
264
 
                }
265
 
              }
266
 
            }
267
 
          }
268
 
        } 
269
 
      }
270
 
      if (is == 1) {
271
 
      // intersection occurs at a box edge - leave lam1 = lam2
272
 
      }
273
 
      else if (is == 2) {
274
 
        lam1 = min(ss[0],ss[1]);
275
 
        lam2 = max(ss[0],ss[1]);
276
 
      }
277
 
      else if (is == 3) {
278
 
      // intersection occurs at a box vertex - leave lam1 = lam2
279
 
      }
280
 
      else {
281
 
        if (is != 0) throw ATC_Error(0,"logic failure in HardyKernel Cell for corner intersection\n");
282
 
      }
283
 
    }
284
 
  }
285
 
 
286
 
  //------------------------------------------------------------------------
287
 
  // constructor
288
 
  ATC_HardyKernelCubicSphere::ATC_HardyKernelCubicSphere
289
 
    (int nparameters, double* parameters): 
290
 
    ATC_HardyKernel(nparameters, parameters) 
291
 
  {
292
 
    for (int k = 0; k < nsd_; k++ ) {
293
 
      if ((bool) periodicity[k]) {
294
 
        if (Rc_ > 0.5*box_length[k]) {
295
 
          throw ATC_Error(0,"Size of localization volume is too large for periodic boundary condition");
296
 
        };
297
 
      };
298
 
    };
299
 
  }
300
 
 
301
 
  // function value
302
 
  double ATC_HardyKernelCubicSphere::value(DENS_VEC& x_atom) 
303
 
  {
304
 
     double r=x_atom.norm();
305
 
     double rn=r/Rc_;
306
 
     if (rn < 1.0) { return 5.0*(1.0-3.0*rn*rn+2.0*rn*rn*rn); }
307
 
     else          { return 0.0; }
308
 
  } 
309
 
  //------------------------------------------------------------------------
310
 
  // constructor
311
 
  ATC_HardyKernelQuarticSphere::ATC_HardyKernelQuarticSphere
312
 
    (int nparameters, double* parameters): 
313
 
    ATC_HardyKernel(nparameters, parameters) 
314
 
  {
315
 
    for (int k = 0; k < nsd_; k++ ) {
316
 
      if ((bool) periodicity[k]) {
317
 
        if (Rc_ > 0.5*box_length[k]) {
318
 
          throw ATC_Error(0,"Size of localization volume is too large for periodic boundary condition");
319
 
        };
320
 
      };
321
 
    };
322
 
  }
323
 
 
324
 
  // function value
325
 
  double ATC_HardyKernelQuarticSphere::value(DENS_VEC& x_atom) 
326
 
  {
327
 
     double r=x_atom.norm();
328
 
     double rn=r/Rc_;
329
 
     if (rn < 1.0) { return 35.0/8.0*pow((1.0-rn*rn),2); }
330
 
     else          { return 0.0; }
331
 
  } 
332
 
  //------------------------------------------------------------------------
333
 
  // constructor
334
 
  ATC_HardyKernelCubicCyl::ATC_HardyKernelCubicCyl
335
 
    (int nparameters, double* parameters): 
336
 
    ATC_HardyKernel(nparameters, parameters)
337
 
  {
338
 
    nsd_ = 2;
339
 
    double Lz = box_length[2];
340
 
    invVol_ = 1.0/(Pi*pow(Rc_,2)*Lz);
341
 
    for (int k = 0; k < nsd_; k++ ) {
342
 
      if ((bool) periodicity[k]) {
343
 
        if (Rc_ > 0.5*box_length[k]) {
344
 
          throw ATC_Error(0,"Size of localization volume is too large for periodic boundary condition");
345
 
        };
346
 
      };
347
 
    };
348
 
  }
349
 
 
350
 
  // function value
351
 
  double ATC_HardyKernelCubicCyl::value(DENS_VEC& x_atom) 
352
 
  {
353
 
     double r=sqrt(pow(x_atom(0),2)+pow(x_atom(1),2));
354
 
     double rn=r/Rc_;
355
 
     if (rn < 1.0) { return 10.0/3.0*(1.0-3.0*rn*rn+2.0*rn*rn*rn); }
356
 
     else          { return 0.0; }
357
 
  } 
358
 
 
359
 
 
360
 
  //------------------------------------------------------------------------
361
 
  // constructor
362
 
  ATC_HardyKernelQuarticCyl::ATC_HardyKernelQuarticCyl
363
 
    (int nparameters, double* parameters): 
364
 
    ATC_HardyKernel(nparameters, parameters)
365
 
  {
366
 
    nsd_ = 2;
367
 
    double Lz = box_length[2];
368
 
    invVol_ = 1.0/(Pi*pow(Rc_,2)*Lz);
369
 
    for (int k = 0; k < nsd_; k++ ) {
370
 
      if ((bool) periodicity[k]) {
371
 
        if (Rc_ > 0.5*box_length[k]) {
372
 
          throw ATC_Error(0,"Size of localization volume is too large for periodic boundary condition");
373
 
        };
374
 
      };
375
 
    };
376
 
  }
377
 
 
378
 
  // function value
379
 
  double ATC_HardyKernelQuarticCyl::value(DENS_VEC& x_atom) 
380
 
  {
381
 
    double r=sqrt(pow(x_atom(0),2)+pow(x_atom(1),2));
382
 
    double rn=r/Rc_;
383
 
    if (rn < 1.0) { return 3.0*pow((1.0-rn*rn),2); }
384
 
    else          { return 0.0; }
385
 
  } 
386
 
};