ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/Electrostatic.cpp
(Generate patch)

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1876 by gezelter, Wed Feb 20 15:39:39 2013 UTC vs.
Revision 1877 by gezelter, Thu Jun 6 15:43:35 2013 UTC

# Line 75 | Line 75 | namespace OpenMD {
75      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
76      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
77      summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;    
78 +    summationMap_["TAYLOR_SHIFTED"]     = esm_TAYLOR_SHIFTED;    
79      summationMap_["REACTION_FIELD"]     = esm_REACTION_FIELD;    
80      summationMap_["EWALD_FULL"]         = esm_EWALD_FULL;        
81      summationMap_["EWALD_PME"]          = esm_EWALD_PME;        
# Line 107 | Line 108 | namespace OpenMD {
108      debyeToCm_ = 3.33564095198e-30;
109      
110      // Default number of points for electrostatic splines
111 <    np_ = 140;
111 >    np_ = 100;
112      
113      // variables to handle different summation methods for long-range
114      // electrostatics:
# Line 129 | Line 130 | namespace OpenMD {
130                   "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
131                   "\t(Input file specified %s .)\n"
132                   "\telectrostaticSummationMethod must be one of: \"hard\",\n"
133 <                 "\t\"shifted_potential\", \"shifted_force\", or \n"
134 <                 "\t\"reaction_field\".\n", myMethod.c_str() );
133 >                 "\t\"shifted_potential\", \"shifted_force\",\n"
134 >                 "\t\"taylor_shifted\", or \"reaction_field\".\n",
135 >                 myMethod.c_str() );
136          painCave.isFatal = 1;
137          simError();
138        }
# Line 234 | Line 236 | namespace OpenMD {
236      
237      RealType r = cutoffRadius_;
238      RealType r2 = r * r;
239 +    RealType ric = 1.0 / r;
240 +    RealType ric2 = ric * ric;
241  
242      if (screeningMethod_ == DAMPED) {      
243        a2 = dampingAlpha_ * dampingAlpha_;
# Line 265 | Line 269 | namespace OpenMD {
269      db0c_4 =          3.0*b2c  - 6.0*r2*b3c     + r2*r2*b4c;
270      db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;
271      
272 +
273      // working variables for the splines:
274      RealType ri, ri2;
275      RealType b0, b1, b2, b3, b4, b5;
276      RealType db0_1, db0_2, db0_3, db0_4, db0_5;
277 <    RealType f0;
278 <    RealType g0, g1, g2, g3, g4;
279 <    RealType h1, h2, h3, h4;
280 <    RealType s2, s3, s4;
281 <    RealType t3, t4;
282 <    RealType u4;
277 >    RealType f, fc, f0;
278 >    RealType g, gc, g0, g1, g2, g3, g4;
279 >    RealType h, hc, h1, h2, h3, h4;
280 >    RealType s, sc, s2, s3, s4;
281 >    RealType t, tc, t3, t4;
282 >    RealType u, uc, u4;
283  
284      // working variables for Taylor expansion:
285      RealType rmRc, rmRc2, rmRc3, rmRc4;
# Line 294 | Line 299 | namespace OpenMD {
299  
300      // Storage vectors for the computed functions    
301      vector<RealType> rv;
302 <    vector<RealType> v01v, v02v;
303 <    vector<RealType> v11v, v12v, v13v;
304 <    vector<RealType> v21v, v22v, v23v, v24v;
305 <    vector<RealType> v31v, v32v, v33v, v34v, v35v;
306 <    vector<RealType> v41v, v42v, v43v, v44v, v45v, v46v;
302 >    vector<RealType> v01v;
303 >    vector<RealType> v11v;
304 >    vector<RealType> v21v, v22v;
305 >    vector<RealType> v31v, v32v;
306 >    vector<RealType> v41v, v42v, v43v;
307  
308 +    /*
309 +    vector<RealType> dv01v;
310 +    vector<RealType> dv11v;
311 +    vector<RealType> dv21v, dv22v;
312 +    vector<RealType> dv31v, dv32v;
313 +    vector<RealType> dv41v, dv42v, dv43v;
314 +    */
315 +
316      for (int i = 1; i < np_ + 1; i++) {
317        r = RealType(i) * dx;
318        rv.push_back(r);
# Line 339 | Line 352 | namespace OpenMD {
352        db0_3 =          3.0*r*b2   - r2*r*b3;
353        db0_4 =          3.0*b2   - 6.0*r2*b3     + r2*r2*b4;
354        db0_5 =                    -15.0*r*b3 + 10.0*r2*r*b4 - r2*r2*r*b5;
355 +
356 +      f = b0;
357 +      fc = b0c;
358 +      f0 = f - fc - rmRc*db0c_1;
359 +
360 +      g = db0_1;        
361 +      gc = db0c_1;
362 +      g0 = g - gc;
363 +      g1 = g0 - rmRc *db0c_2;
364 +      g2 = g1 - rmRc2*db0c_3;
365 +      g3 = g2 - rmRc3*db0c_4;
366 +      g4 = g3 - rmRc4*db0c_5;
367 +
368 +      h = db0_2;      
369 +      hc = db0c_2;
370 +      h1 = h - hc;
371 +      h2 = h1 - rmRc *db0c_3;
372 +      h3 = h2 - rmRc2*db0c_4;
373 +      h4 = h3 - rmRc3*db0c_5;
374 +
375 +      s = db0_3;      
376 +      sc = db0c_3;
377 +      s2 = s - sc;
378 +      s3 = s2 - rmRc *db0c_4;
379 +      s4 = s3 - rmRc2*db0c_5;
380  
381 +      t = db0_4;      
382 +      tc = db0c_4;
383 +      t3 = t - tc;
384 +      t4 = t3 - rmRc *db0c_5;
385 +      
386 +      u = db0_5;        
387 +      uc = db0c_5;
388 +      u4 = u - uc;
389  
390 +      // in what follows below, the various v functions are used for
391 +      // potentials and torques, while the w functions show up in the
392 +      // forces.
393 +
394        switch (summationMethod_) {
395        case esm_SHIFTED_FORCE:
396 <        f0 = b0 - b0c - rmRc*db0c_1;
396 >                
397 >        v01 = f - fc - rmRc*gc;
398 >        v11 = g - gc - rmRc*hc;
399 >        v21 = g*ri - gc*ric - rmRc*(hc - gc*ric)*ric;
400 >        v22 = h - g*ri - (hc - gc*ric) - rmRc*(sc - (hc - gc*ric)*ric);
401 >        v31 = (h-g*ri)*ri - (hc-g*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
402 >        v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric)
403 >          - rmRc*(tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
404 >        v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2
405 >          - rmRc*(sc - 3.0*(hc-gc*ric)*ric)*ric2;
406 >        v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric
407 >          - rmRc*(tc - (4.0*sc - 9.0*(hc - gc*ric)*ric)*ric)*ric;
408          
409 <        g0 = db0_1 - db0c_1;
410 <        g1 = g0 - rmRc *db0c_2;
411 <        g2 = g1 - rmRc2*db0c_3;
412 <        g3 = g2 - rmRc3*db0c_4;
413 <        g4 = g3 - rmRc4*db0c_5;
409 >        v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
410 >          - (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric)
411 >          - rmRc*(uc-3.0*(2.0*tc - (7.0*sc - 15.0*(hc - gc*ric)*ric)*ric)*ric);
412 >
413 >        dv01 = g - gc;
414 >        dv11 = h - hc;
415 >        dv21 = (h - g*ri)*ri - (hc - gc*ric)*ric;
416 >        dv22 = (s - (h - g*ri)*ri) - (sc - (hc - gc*ric)*ric);        
417 >        dv31 = (s - 2.0*(h-g*ri)*ri)*ri - (sc - 2.0*(hc-gc*ric)*ric)*ric;
418 >        dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri)
419 >          - (tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
420 >        dv41 = (s - 3.0*(h - g*ri)*ri)*ri2 - (sc - 3.0*(hc - gc*ric)*ric)*ric2;
421 >        dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri
422 >          - (tc - (4.0*sc - 9.0*(hc-gc*ric)*ric)*ric)*ric;
423 >        dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri)
424 >          - (uc - 3.0*(2.0*tc - (7.0*sc - 15.0*(hc - gc*ric)*ric)*ric)*ric);
425          
426 <        h1 = db0_2 - db0c_2;
427 <        h2 = h1 - rmRc *db0c_3;
428 <        h3 = h2 - rmRc2*db0c_4;
357 <        h4 = h3 - rmRc3*db0c_5;
426 >        break;
427 >
428 >      case esm_TAYLOR_SHIFTED:
429          
430 <        s2 = db0_3 - db0c_3;
431 <        s3 = s2 - rmRc *db0c_4;
432 <        s4 = s3 - rmRc2*db0c_5;
433 <        
434 <        t3 = db0_4 - db0c_4;
435 <        t4 = t3 - rmRc *db0c_5;
436 <        
437 <        u4 = db0_5 - db0c_5;
430 >        v01 = f0;
431 >        v11 = g1;
432 >        v21 = g2 * ri;
433 >        v22 = h2 - v21;
434 >        v31 = (h3 - g3 * ri) * ri;
435 >        v32 = s3 - 3.0*v31;
436 >        v41 = (h4 - g4 * ri) * ri2;
437 >        v42 = s4 * ri - 3.0*v41;
438 >        v43 = t4 - 6.0*v42 - 3.0*v41;
439 >
440 >        dv01 = g0;
441 >        dv11 = h1;
442 >        dv21 = (h2 - g2*ri)*ri;
443 >        dv22 = (s2 - (h2 - g2*ri)*ri);
444 >        dv31 = (s3 - 2.0*(h3-g3*ri)*ri)*ri;
445 >        dv32 = (t3 - 3.0*(s3-2.0*(h3-g3*ri)*ri)*ri);
446 >        dv41 = (s4 - 3.0*(h4 - g4*ri)*ri)*ri2;
447 >        dv42 = (t4 - (4.0*s4 - 9.0*(h4-g4*ri)*ri)*ri)*ri;
448 >        dv43 = (u4 - 3.0*(2.0*t4 - (7.0*s4 - 15.0*(h4 - g4*ri)*ri)*ri)*ri);
449 >
450          break;
451  
452        case esm_SHIFTED_POTENTIAL:
370 –        f0 = b0 - b0c;
371 –        
372 –        g0 = db0_1;
373 –        g1 = db0_1 - db0c_1;
374 –        g2 = g1 - rmRc *db0c_2;
375 –        g3 = g2 - rmRc2*db0c_3;
376 –        g4 = g3 - rmRc3*db0c_4;
453  
454 <        h1 = db0_2;
455 <        h2 = db0_2 - db0c_2;
456 <        h3 = h2 - rmRc *db0c_3;
457 <        h4 = h3 - rmRc2*db0c_4;
458 <        
459 <        s2 = db0_3;
460 <        s3 = db0_3 - db0c_3;
461 <        s4 = s3 - rmRc *db0c_4;
454 >        v01 = f - fc;
455 >        v11 = g - gc;
456 >        v21 = g*ri - gc*ric;
457 >        v22 = h - g*ri - (hc - gc*ric);
458 >        v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
459 >        v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
460 >        v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
461 >        v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;        
462 >        v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
463 >          - (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric);
464  
465 <        t3 = db0_4;
466 <        t4 = db0_4 - db0c_4;
467 <        
468 <        u4 = db0_5;
465 >        dv01 = g;
466 >        dv11 = h;
467 >        dv21 = (h - g*ri)*ri;
468 >        dv22 = (s - (h - g*ri)*ri);
469 >        dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
470 >        dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
471 >        dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
472 >        dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
473 >        dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
474 >
475          break;
476  
477        case esm_SWITCHING_FUNCTION:
478        case esm_HARD:
479 <        f0 = b0;
480 <        
481 <        g0 = db0_1;
482 <        g1 = g0;
483 <        g2 = g1;
484 <        g3 = g2;
485 <        g4 = g3;
486 <        
487 <        h1 = db0_2;
488 <        h2 = h1;
489 <        h3 = h2;
490 <        h4 = h3;
491 <        
492 <        s2 = db0_3;
493 <        s3 = s2;
494 <        s4 = s3;
495 <        
496 <        t3 = db0_4;
497 <        t4 = t3;
498 <        
499 <        u4 = db0_5;
479 >
480 >        v01 = f;
481 >        v11 = g;
482 >        v21 = g*ri;
483 >        v22 = h - g*ri;
484 >        v31 = (h-g*ri)*ri;
485 >        v32 = (s - 3.0*(h-g*ri)*ri);
486 >        v41 = (h - g*ri)*ri2;
487 >        v42 = (s-3.0*(h-g*ri)*ri)*ri;        
488 >        v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri);
489 >
490 >        dv01 = g;
491 >        dv11 = h;
492 >        dv21 = (h - g*ri)*ri;
493 >        dv22 = (s - (h - g*ri)*ri);
494 >        dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
495 >        dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
496 >        dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
497 >        dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
498 >        dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
499 >
500          break;
501  
502        case esm_REACTION_FIELD:
503 <
503 >        
504          // following DL_POLY's lead for shifting the image charge potential:
505 <        f0 = b0  + preRF_ * r2
506 <          - (b0c + preRF_ * cutoffRadius_ * cutoffRadius_);
505 >        f = b0 + preRF_ * r2;
506 >        fc = b0c + preRF_ * cutoffRadius_ * cutoffRadius_;
507  
508 <        g0 = db0_1 + preRF_ * 2.0 * r;
509 <        g1 = g0;
426 <        g2 = g1;
427 <        g3 = g2;
428 <        g4 = g3;
429 <
430 <        h1 = db0_2 + preRF_ * 2.0;
431 <        h2 = h1;
432 <        h3 = h2;
433 <        h4 = h3;
508 >        g = db0_1 + preRF_ * 2.0 * r;        
509 >        gc = db0c_1 + preRF_ * 2.0 * cutoffRadius_;
510  
511 <        s2 = db0_3;
512 <        s3 = s2;
513 <        s4 = s3;
514 <        
515 <        t3 = db0_4;
516 <        t4 = t3;
517 <        
518 <        u4 = db0_5;        
511 >        h = db0_2 + preRF_ * 2.0;
512 >        hc = db0c_2 + preRF_ * 2.0;
513 >
514 >        v01 = f - fc;
515 >        v11 = g - gc;
516 >        v21 = g*ri - gc*ric;
517 >        v22 = h - g*ri - (hc - gc*ric);
518 >        v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
519 >        v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
520 >        v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
521 >        v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;        
522 >        v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
523 >          - (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric);
524 >
525 >        dv01 = g;
526 >        dv11 = h;
527 >        dv21 = (h - g*ri)*ri;
528 >        dv22 = (s - (h - g*ri)*ri);
529 >        dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
530 >        dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
531 >        dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
532 >        dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
533 >        dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
534 >
535          break;
536                  
537        case esm_EWALD_FULL:
# Line 461 | Line 553 | namespace OpenMD {
553          break;      
554        }
555  
464 –      v01 = f0;
465 –      v02 = g0;
466 –
467 –      v11 = g1;
468 –      v12 = g1 * ri;
469 –      v13 = h1 - v12;
470 –
471 –      v21 = g2 * ri;
472 –      v22 = h2 - v21;
473 –      v23 = v22 * ri;
474 –      v24 = s2 - 3.0*v23;        
475 –
476 –      v31 = (h3 - g3 * ri) * ri;
477 –      v32 = s3 - 3.0*v31;
478 –      v33 = v31 * ri;
479 –      v34 = v32 * ri;
480 –      v35 = t3 - 6.0*v34 - 3.0*v33;
481 –
482 –      v41 = (h4 - g4 * ri) * ri2;
483 –      v42 = s4 * ri - 3.0*v41;
484 –      v43 = t4 - 6.0*v42 - 3.0*v41;
485 –      v44 = v42 * ri;
486 –      v45 = v43 * ri;
487 –      v46 = u4 - 10.0*v45 - 15.0*v44;
488 –
556        // Add these computed values to the storage vectors for spline creation:
557        v01v.push_back(v01);
491 –      v02v.push_back(v02);
492 –
558        v11v.push_back(v11);
494 –      v12v.push_back(v12);
495 –      v13v.push_back(v13);
496 –
559        v21v.push_back(v21);
560        v22v.push_back(v22);
499 –      v23v.push_back(v23);
500 –      v24v.push_back(v24);
501 –
561        v31v.push_back(v31);
562 <      v32v.push_back(v32);
504 <      v33v.push_back(v33);
505 <      v34v.push_back(v34);
506 <      v35v.push_back(v35);
507 <      
562 >      v32v.push_back(v32);      
563        v41v.push_back(v41);
564        v42v.push_back(v42);
565        v43v.push_back(v43);
566 <      v44v.push_back(v44);
567 <      v45v.push_back(v45);
568 <      v46v.push_back(v46);
566 >      /*
567 >      dv01v.push_back(dv01);
568 >      dv11v.push_back(dv11);
569 >      dv21v.push_back(dv21);
570 >      dv22v.push_back(dv22);
571 >      dv31v.push_back(dv31);
572 >      dv32v.push_back(dv32);      
573 >      dv41v.push_back(dv41);
574 >      dv42v.push_back(dv42);
575 >      dv43v.push_back(dv43);
576 >      */
577      }
578  
579      // construct the spline structures and fill them with the values we've
# Line 518 | Line 581 | namespace OpenMD {
581  
582      v01s = new CubicSpline();
583      v01s->addPoints(rv, v01v);
521 –    v02s = new CubicSpline();
522 –    v02s->addPoints(rv, v02v);
523 –
584      v11s = new CubicSpline();
585      v11s->addPoints(rv, v11v);
526 –    v12s = new CubicSpline();
527 –    v12s->addPoints(rv, v12v);
528 –    v13s = new CubicSpline();
529 –    v13s->addPoints(rv, v13v);
530 –
586      v21s = new CubicSpline();
587      v21s->addPoints(rv, v21v);
588      v22s = new CubicSpline();
589      v22s->addPoints(rv, v22v);
535 –    v23s = new CubicSpline();
536 –    v23s->addPoints(rv, v23v);
537 –    v24s = new CubicSpline();
538 –    v24s->addPoints(rv, v24v);
539 –
590      v31s = new CubicSpline();
591      v31s->addPoints(rv, v31v);
592      v32s = new CubicSpline();
593      v32s->addPoints(rv, v32v);
544 –    v33s = new CubicSpline();
545 –    v33s->addPoints(rv, v33v);
546 –    v34s = new CubicSpline();
547 –    v34s->addPoints(rv, v34v);
548 –    v35s = new CubicSpline();
549 –    v35s->addPoints(rv, v35v);
550 –
594      v41s = new CubicSpline();
595      v41s->addPoints(rv, v41v);
596      v42s = new CubicSpline();
597      v42s->addPoints(rv, v42v);
598      v43s = new CubicSpline();
599      v43s->addPoints(rv, v43v);
557 –    v44s = new CubicSpline();
558 –    v44s->addPoints(rv, v44v);
559 –    v45s = new CubicSpline();
560 –    v45s->addPoints(rv, v45v);
561 –    v46s = new CubicSpline();
562 –    v46s->addPoints(rv, v46v);
600  
601 +    /*
602 +    dv01s = new CubicSpline();
603 +    dv01s->addPoints(rv, dv01v);
604 +    dv11s = new CubicSpline();
605 +    dv11s->addPoints(rv, dv11v);
606 +    dv21s = new CubicSpline();
607 +    dv21s->addPoints(rv, dv21v);
608 +    dv22s = new CubicSpline();
609 +    dv22s->addPoints(rv, dv22v);
610 +    dv31s = new CubicSpline();
611 +    dv31s->addPoints(rv, dv31v);
612 +    dv32s = new CubicSpline();
613 +    dv32s->addPoints(rv, dv32v);
614 +    dv41s = new CubicSpline();
615 +    dv41s->addPoints(rv, dv41v);
616 +    dv42s = new CubicSpline();
617 +    dv42s->addPoints(rv, dv42v);
618 +    dv43s = new CubicSpline();
619 +    dv43s->addPoints(rv, dv43v);
620 +    */
621 +
622      haveElectroSplines_ = true;
623  
624      initialized_ = true;
# Line 753 | Line 811 | namespace OpenMD {
811      
812      // needed for fields (and forces):
813      if (a_is_Charge || b_is_Charge) {
814 <      v02 = v02s->getValueAt( *(idat.rij) );
814 >      v01s->getValueAndDerivativeAt( *(idat.rij), v01, dv01);
815      }
816      if (a_is_Dipole || b_is_Dipole) {
817 <      v12 = v12s->getValueAt( *(idat.rij) );
818 <      v13 = v13s->getValueAt( *(idat.rij) );
817 >      v11s->getValueAndDerivativeAt( *(idat.rij), v11, dv11);
818 >      v11or = ri * v11;
819      }
820 <    if (a_is_Quadrupole || b_is_Quadrupole) {
821 <      v23 = v23s->getValueAt( *(idat.rij) );
822 <      v24 = v24s->getValueAt( *(idat.rij) );
823 <    }
820 >    if (a_is_Quadrupole || b_is_Quadrupole ||  (a_is_Dipole && b_is_Dipole)) {
821 >      v21s->getValueAndDerivativeAt( *(idat.rij), v21, dv21);
822 >      v22s->getValueAndDerivativeAt( *(idat.rij), v22, dv22);
823 >      v22or = ri * v22;
824 >    }      
825  
826 <    // needed for potentials (and torques):
768 <    if (a_is_Charge && b_is_Charge) {
769 <      v01 = v01s->getValueAt( *(idat.rij) );
770 <    }
771 <    if ((a_is_Charge && b_is_Dipole) || (b_is_Charge && a_is_Dipole)) {
772 <      v11 = v11s->getValueAt( *(idat.rij) );
773 <    }
774 <    if ((a_is_Charge && b_is_Quadrupole) || (b_is_Charge && a_is_Quadrupole)) {
775 <      v21 = v21s->getValueAt( *(idat.rij) );
776 <      v22 = v22s->getValueAt( *(idat.rij) );
777 <    } else if (a_is_Dipole && b_is_Dipole) {
778 <      v21 = v21s->getValueAt( *(idat.rij) );
779 <      v22 = v22s->getValueAt( *(idat.rij) );
780 <      v23 = v23s->getValueAt( *(idat.rij) );
781 <      v24 = v24s->getValueAt( *(idat.rij) );
782 <    }
826 >    // needed for potentials (and forces and torques):
827      if ((a_is_Dipole && b_is_Quadrupole) ||
828          (b_is_Dipole && a_is_Quadrupole)) {
829 <      v31 = v31s->getValueAt( *(idat.rij) );
830 <      v32 = v32s->getValueAt( *(idat.rij) );
831 <      v33 = v33s->getValueAt( *(idat.rij) );
832 <      v34 = v34s->getValueAt( *(idat.rij) );
789 <      v35 = v35s->getValueAt( *(idat.rij) );
829 >      v31s->getValueAndDerivativeAt( *(idat.rij), v31, dv31);
830 >      v32s->getValueAndDerivativeAt( *(idat.rij), v32, dv32);
831 >      v31or = v31 * ri;
832 >      v32or = v32 * ri;
833      }
834      if (a_is_Quadrupole && b_is_Quadrupole) {
835 <      v41 = v41s->getValueAt( *(idat.rij) );
836 <      v42 = v42s->getValueAt( *(idat.rij) );
837 <      v43 = v43s->getValueAt( *(idat.rij) );
838 <      v44 = v44s->getValueAt( *(idat.rij) );
839 <      v45 = v45s->getValueAt( *(idat.rij) );
797 <      v46 = v46s->getValueAt( *(idat.rij) );
835 >      v41s->getValueAndDerivativeAt( *(idat.rij), v41, dv41);
836 >      v42s->getValueAndDerivativeAt( *(idat.rij), v42, dv42);
837 >      v43s->getValueAndDerivativeAt( *(idat.rij), v43, dv43);
838 >      v42or = v42 * ri;
839 >      v43or = v43 * ri;
840      }
841  
842      // calculate the single-site contributions (fields, etc).
# Line 810 | Line 852 | namespace OpenMD {
852          *(idat.skippedCharge2) += C_a;
853        } else {
854          // only do the field if we're not excluded:
855 <        Eb -= C_a *  pre11_ * v02 * rhat;
855 >        Eb -= C_a *  pre11_ * dv01 * rhat;
856        }
857      }
858      
# Line 819 | Line 861 | namespace OpenMD {
861        rdDa = dot(rhat, D_a);
862        rxDa = cross(rhat, D_a);
863        if (!idat.excluded)
864 <        Eb -=  pre12_ * (v13 * rdDa * rhat + v12 * D_a);
864 >        Eb -=  pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
865      }
866      
867      if (a_is_Quadrupole) {
# Line 830 | Line 872 | namespace OpenMD {
872        rdQar = dot(rhat, Qar);
873        rxQar = cross(rhat, Qar);
874        if (!idat.excluded)
875 <        Eb -= pre14_ * ((trQa * rhat + 2.0 * Qar) * v23 + rdQar * rhat * v24);
875 >        Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or
876 >                        + rdQar * rhat * (dv22 - 2.0*v22or));
877      }
878      
879      if (b_is_Charge) {
# Line 843 | Line 886 | namespace OpenMD {
886          *(idat.skippedCharge1) += C_b;
887        } else {
888          // only do the field if we're not excluded:
889 <        Ea += C_b *  pre11_ * v02 * rhat;
889 >        Ea += C_b *  pre11_ * dv01 * rhat;
890        }
891      }
892      
# Line 852 | Line 895 | namespace OpenMD {
895        rdDb = dot(rhat, D_b);
896        rxDb = cross(rhat, D_b);
897        if (!idat.excluded)
898 <        Ea += pre12_ * (v13 * rdDb * rhat + v12 * D_b);
898 >        Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b);
899      }
900      
901      if (b_is_Quadrupole) {
# Line 863 | Line 906 | namespace OpenMD {
906        rdQbr = dot(rhat, Qbr);
907        rxQbr = cross(rhat, Qbr);
908        if (!idat.excluded)
909 <        Ea += pre14_ * ((trQb * rhat + 2.0 * Qbr) * v23 + rdQbr * rhat * v24);
909 >        Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
910 >                        + rdQbr * rhat * (dv22 - 2.0*v22or));
911      }
912      
913      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
# Line 873 | Line 917 | namespace OpenMD {
917      if (a_is_Charge) {    
918        
919        if (b_is_Charge) {
920 <        pref =  pre11_ * *(idat.electroMult);          
920 >        pref =  pre11_ * *(idat.electroMult);      
921          U  += C_a * C_b * pref * v01;
922 <        F  += C_a * C_b * pref * v02 * rhat;
922 >        F  += C_a * C_b * pref * dv01 * rhat;
923          
924          // If this is an excluded pair, there are still indirect
925          // interactions via the reaction field we must worry about:
# Line 906 | Line 950 | namespace OpenMD {
950        if (b_is_Dipole) {
951          pref =  pre12_ * *(idat.electroMult);        
952          U  += C_a * pref * v11 * rdDb;
953 <        F  += C_a * pref * (v13 * rdDb * rhat + v12 * D_b);
953 >        F  += C_a * pref * ((dv11 - v11or) * rdDb * rhat + v11or * D_b);
954          Tb += C_a * pref * v11 * rxDb;
955  
956          if (a_is_Fluctuating) dUdCa += pref * v11 * rdDb;
# Line 926 | Line 970 | namespace OpenMD {
970        if (b_is_Quadrupole) {
971          pref = pre14_ * *(idat.electroMult);
972          U  +=  C_a * pref * (v21 * trQb + v22 * rdQbr);
973 <        F  +=  C_a * pref * (trQb * rhat + 2.0 * Qbr) * v23;
974 <        F  +=  C_a * pref * rdQbr * rhat * v24;
973 >        F  +=  C_a * pref * (trQb * dv21 * rhat + 2.0 * Qbr * v22or);
974 >        F  +=  C_a * pref * rdQbr * rhat * (dv22 - 2.0*v22or);
975          Tb +=  C_a * pref * 2.0 * rxQbr * v22;
976  
977          if (a_is_Fluctuating) dUdCa += pref * (v21 * trQb + v22 * rdQbr);
# Line 940 | Line 984 | namespace OpenMD {
984          pref = pre12_ * *(idat.electroMult);
985  
986          U  -= C_b * pref * v11 * rdDa;
987 <        F  -= C_b * pref * (v13 * rdDa * rhat + v12 * D_a);
987 >        F  -= C_b * pref * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
988          Ta -= C_b * pref * v11 * rxDa;
989  
990          if (b_is_Fluctuating) dUdCb -= pref * v11 * rdDa;
# Line 962 | Line 1006 | namespace OpenMD {
1006          DaxDb = cross(D_a, D_b);
1007  
1008          U  -= pref * (DadDb * v21 + rdDa * rdDb * v22);
1009 <        F  -= pref * (DadDb * rhat + rdDb * D_a + rdDa * D_b)*v23;
1010 <        F  -= pref * (rdDa * rdDb) * v24 * rhat;
1009 >        F  -= pref * (dv21 * DadDb * rhat + v22or * (rdDb * D_a + rdDa * D_b));
1010 >        F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
1011          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
1012          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
1013  
# Line 976 | Line 1020 | namespace OpenMD {
1020            indirect_Ta  += rfContrib * DaxDb;
1021            indirect_Tb  -= rfContrib * DaxDb;
1022          }
979 –
1023        }
1024  
1025        if (b_is_Quadrupole) {
# Line 986 | Line 1029 | namespace OpenMD {
1029          DaxQbr = cross(D_a, Qbr);
1030  
1031          U  -= pref * ((trQb*rdDa + 2.0*DadQbr)*v31 + rdDa*rdQbr*v32);
1032 <        F  -= pref * (trQb*D_a + 2.0*DadQb) * v33;
1033 <        F  -= pref * (trQb*rdDa*rhat + 2.0*DadQbr*rhat + D_a*rdQbr
1034 <                      + 2.0*rdDa*rQb)*v34;
1035 <        F  -= pref * (rdDa * rdQbr * rhat * v35);
1032 >        F  -= pref * (trQb*D_a + 2.0*DadQb) * v31or;
1033 >        F  -= pref * (trQb*rdDa + 2.0*DadQbr) * (dv31-v31or) * rhat;
1034 >        F  -= pref * (D_a*rdQbr + 2.0*rdDa*rQb) * v32or;
1035 >        F  -= pref * (rdDa * rdQbr * rhat * (dv32-3.0*v32or));
1036          Ta += pref * ((-trQb*rxDa + 2.0 * DaxQbr)*v31 - rxDa*rdQbr*v32);
1037          Tb += pref * ((2.0*cross(DadQb, rhat) - 2.0*DaxQbr)*v31
1038                        - 2.0*rdDa*rxQbr*v32);
# Line 1000 | Line 1043 | namespace OpenMD {
1043        if (b_is_Charge) {
1044          pref = pre14_ * *(idat.electroMult);
1045          U  += C_b * pref * (v21 * trQa + v22 * rdQar);
1046 <        F  += C_b * pref * (trQa * rhat + 2.0 * Qar) * v23;
1047 <        F  += C_b * pref * rdQar * rhat * v24;
1046 >        F  += C_b * pref * (trQa * rhat * dv21 + 2.0 * Qar * v22or);
1047 >        F  += C_b * pref * rdQar * rhat * (dv22 - 2.0*v22or);
1048          Ta += C_b * pref * 2.0 * rxQar * v22;
1049  
1050          if (b_is_Fluctuating) dUdCb += pref * (v21 * trQa + v22 * rdQar);
# Line 1013 | Line 1056 | namespace OpenMD {
1056          DbxQar = cross(D_b, Qar);
1057  
1058          U  += pref * ((trQa*rdDb + 2.0*DbdQar)*v31 + rdDb*rdQar*v32);
1059 <        F  += pref * (trQa*D_b + 2.0*DbdQa) * v33;
1060 <        F  += pref * (trQa*rdDb*rhat + 2.0*DbdQar*rhat + D_b*rdQar
1061 <                      + 2.0*rdDb*rQa)*v34;
1062 <        F  += pref * (rdDb * rdQar * rhat * v35);
1059 >        F  += pref * (trQa*D_b + 2.0*DbdQa) * v31or;
1060 >        F  += pref * (trQa*rdDb + 2.0*DbdQar) * (dv31-v31or) * rhat;
1061 >        F  += pref * (D_b*rdQar + 2.0*rdDb*rQa) * v32or;
1062 >        F  += pref * (rdDb * rdQar * rhat * (dv32-3.0*v32or));
1063          Ta += pref * ((-2.0*cross(DbdQa, rhat) + 2.0*DbxQar)*v31
1064                        + 2.0*rdDb*rxQar*v32);
1065          Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
# Line 1035 | Line 1078 | namespace OpenMD {
1078          U  += pref * (trQa * rdQbr + trQb * rdQar  + 4.0 * rQaQbr) * v42;
1079          U  += pref * (rdQar * rdQbr) * v43;
1080  
1081 <        F  += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*v44;
1082 <        F  += pref * rhat * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr)*v45;
1083 <        F  += pref * rhat * (rdQar * rdQbr)*v46;
1081 >        F  += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*dv41;
1082 >        F  += pref*rhat*(trQa*rdQbr + trQb*rdQar + 4.0*rQaQbr)*(dv42-2.0*v42or);
1083 >        F  += pref * rhat * (rdQar * rdQbr)*(dv43 - 4.0*v43or);
1084  
1085 <        F  += pref * 2.0 * (trQb*rQa + trQa*rQb)*v44;
1086 <        F  += pref * 4.0 * (rQaQb + QaQbr)*v44;
1087 <        F  += pref * 2.0 * (rQa*rdQbr + rdQar*rQb)*v45;
1085 >        F  += pref * 2.0 * (trQb*rQa + trQa*rQb) * v42or;
1086 >        F  += pref * 4.0 * (rQaQb + QaQbr) * v42or;
1087 >        F  += pref * 2.0 * (rQa*rdQbr + rdQar*rQb) * v43or;
1088  
1089          Ta += pref * (- 4.0 * QaxQb * v41);
1090          Ta += pref * (- 2.0 * trQb * cross(rQa, rhat)

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines