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

Comparing trunk/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 1938 by gezelter, Thu Oct 31 15:32:17 2013 UTC

# Line 40 | Line 40
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43 + #ifdef IS_MPI
44 + #include <mpi.h>
45 + #endif
46 +
47   #include <stdio.h>
48   #include <string.h>
49  
50   #include <cmath>
51 + #include <numeric>
52   #include "nonbonded/Electrostatic.hpp"
53   #include "utils/simError.h"
54   #include "types/NonBondedInteractionType.hpp"
# Line 55 | Line 60
60   #include "utils/PhysicalConstants.hpp"
61   #include "math/erfc.hpp"
62   #include "math/SquareMatrix.hpp"
63 + #include "primitives/Molecule.hpp"
64  
65   namespace OpenMD {
66    
# Line 191 | Line 197 | namespace OpenMD {
197        simError();
198      }
199            
200 <    if (screeningMethod_ == DAMPED) {      
200 >    if (screeningMethod_ == DAMPED || summationMethod_ == esm_EWALD_FULL) {
201        if (!simParams_->haveDampingAlpha()) {
202          // first set a cutoff dependent alpha value
203          // we assume alpha depends linearly with rcut from 0 to 20.5 ang
204          dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
205 <        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
200 <        
205 >        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;        
206          // throw warning
207          sprintf( painCave.errMsg,
208                   "Electrostatic::initialize: dampingAlpha was not specified in the\n"
# Line 213 | Line 218 | namespace OpenMD {
218        haveDampingAlpha_ = true;
219      }
220  
221 <    // find all of the Electrostatic atom Types:
222 <    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
223 <    ForceField::AtomTypeContainer::MapTypeIterator i;
224 <    AtomType* at;
221 >
222 >    Etypes.clear();
223 >    Etids.clear();
224 >    FQtypes.clear();
225 >    FQtids.clear();
226 >    ElectrostaticMap.clear();
227 >    Jij.clear();
228 >    nElectro_ = 0;
229 >    nFlucq_ = 0;
230 >
231 >    Etids.resize( forceField_->getNAtomType(), -1);
232 >    FQtids.resize( forceField_->getNAtomType(), -1);
233 >
234 >    set<AtomType*>::iterator at;
235 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {    
236 >      if ((*at)->isElectrostatic()) nElectro_++;
237 >      if ((*at)->isFluctuatingCharge()) nFlucq_++;
238 >    }
239      
240 <    for (at = atomTypes->beginType(i); at != NULL;
241 <         at = atomTypes->nextType(i)) {
242 <      
243 <      if (at->isElectrostatic())
225 <        addType(at);
240 >    Jij.resize(nFlucq_);
241 >
242 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
243 >      if ((*at)->isElectrostatic()) addType(*at);
244      }  
245      
246      if (summationMethod_ == esm_REACTION_FIELD) {
# Line 250 | Line 268 | namespace OpenMD {
268        b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
269        b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
270        b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
271 <      selfMult_ = b0c + a2 * invArootPi;
271 >      // Half the Smith self piece:
272 >      selfMult1_ = - a2 * invArootPi;
273 >      selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0;
274 >      selfMult4_ = - 4.0 * a2 * a2 * a2 * invArootPi / 5.0;
275      } else {
276        a2 = 0.0;
277        b0c = 1.0 / r;
# Line 259 | Line 280 | namespace OpenMD {
280        b3c = (5.0 * b2c) / r2;
281        b4c = (7.0 * b3c) / r2;
282        b5c = (9.0 * b4c) / r2;
283 <      selfMult_ = b0c;
283 >      selfMult1_ = 0.0;
284 >      selfMult2_ = 0.0;
285 >      selfMult4_ = 0.0;
286      }
287  
288      // higher derivatives of B_0 at the cutoff radius:
# Line 267 | Line 290 | namespace OpenMD {
290      db0c_2 =     -b1c + r2 * b2c;
291      db0c_3 =          3.0*r*b2c  - r2*r*b3c;
292      db0c_4 =          3.0*b2c  - 6.0*r2*b3c     + r2*r2*b4c;
293 <    db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;
271 <    
293 >    db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;  
294  
295 +    if (summationMethod_ != esm_EWALD_FULL) {
296 +      selfMult1_ -= b0c;
297 +      selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) /  3.0;
298 +      selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0;
299 +    }
300 +
301      // working variables for the splines:
302      RealType ri, ri2;
303      RealType b0, b1, b2, b3, b4, b5;
# Line 305 | Line 333 | namespace OpenMD {
333      vector<RealType> v31v, v32v;
334      vector<RealType> v41v, v42v, v43v;
335  
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 –
336      for (int i = 1; i < np_ + 1; i++) {
337        r = RealType(i) * dx;
338        rv.push_back(r);
# Line 398 | Line 418 | namespace OpenMD {
418          v11 = g - gc - rmRc*hc;
419          v21 = g*ri - gc*ric - rmRc*(hc - gc*ric)*ric;
420          v22 = h - g*ri - (hc - gc*ric) - rmRc*(sc - (hc - gc*ric)*ric);
421 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
421 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
422          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric)
423            - rmRc*(tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
424          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2
# Line 455 | Line 475 | namespace OpenMD {
475          v11 = g - gc;
476          v21 = g*ri - gc*ric;
477          v22 = h - g*ri - (hc - gc*ric);
478 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
478 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric;
479          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
480          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
481          v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;        
# Line 476 | Line 496 | namespace OpenMD {
496  
497        case esm_SWITCHING_FUNCTION:
498        case esm_HARD:
499 +      case esm_EWALD_FULL:
500  
501          v01 = f;
502          v11 = g;
# Line 515 | Line 536 | namespace OpenMD {
536          v11 = g - gc;
537          v21 = g*ri - gc*ric;
538          v22 = h - g*ri - (hc - gc*ric);
539 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
539 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric;
540          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
541          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
542          v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;        
# Line 534 | Line 555 | namespace OpenMD {
555  
556          break;
557                  
537 –      case esm_EWALD_FULL:
558        case esm_EWALD_PME:
559        case esm_EWALD_SPME:
560        default :
# Line 563 | Line 583 | namespace OpenMD {
583        v41v.push_back(v41);
584        v42v.push_back(v42);
585        v43v.push_back(v43);
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 –      */
586      }
587  
588      // construct the spline structures and fill them with the values we've
# Line 598 | Line 607 | namespace OpenMD {
607      v43s = new CubicSpline();
608      v43s->addPoints(rv, v43v);
609  
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 –
610      haveElectroSplines_ = true;
611  
612      initialized_ = true;
613    }
614        
615    void Electrostatic::addType(AtomType* atomType){
616 <
616 >    
617      ElectrostaticAtomData electrostaticAtomData;
618      electrostaticAtomData.is_Charge = false;
619      electrostaticAtomData.is_Dipole = false;
# Line 661 | Line 649 | namespace OpenMD {
649        electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
650      }
651  
652 <    pair<map<int,AtomType*>::iterator,bool> ret;    
653 <    ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
654 <                                                        atomType) );
652 >    int atid = atomType->getIdent();
653 >    int etid = Etypes.size();
654 >    int fqtid = FQtypes.size();
655 >
656 >    pair<set<int>::iterator,bool> ret;    
657 >    ret = Etypes.insert( atid );
658      if (ret.second == false) {
659        sprintf( painCave.errMsg,
660                 "Electrostatic already had a previous entry with ident %d\n",
661 <               atomType->getIdent() );
661 >               atid);
662        painCave.severity = OPENMD_INFO;
663        painCave.isFatal = 0;
664        simError();        
665      }
666      
667 <    ElectrostaticMap[atomType] = electrostaticAtomData;  
667 >    Etids[ atid ] = etid;
668 >    ElectrostaticMap.push_back(electrostaticAtomData);
669  
670 <    // Now, iterate over all known types and add to the mixing map:
671 <    
672 <    map<AtomType*, ElectrostaticAtomData>::iterator it;
673 <    for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
674 <      AtomType* atype2 = (*it).first;
675 <      ElectrostaticAtomData eaData2 = (*it).second;
676 <      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
677 <        
670 >    if (electrostaticAtomData.is_Fluctuating) {
671 >      ret = FQtypes.insert( atid );
672 >      if (ret.second == false) {
673 >        sprintf( painCave.errMsg,
674 >                 "Electrostatic already had a previous fluctuating charge entry with ident %d\n",
675 >                 atid );
676 >        painCave.severity = OPENMD_INFO;
677 >        painCave.isFatal = 0;
678 >        simError();        
679 >      }
680 >      FQtids[atid] = fqtid;
681 >      Jij[fqtid].resize(nFlucq_);
682 >
683 >      // Now, iterate over all known fluctuating and add to the
684 >      // coulomb integral map:
685 >      
686 >      std::set<int>::iterator it;
687 >      for( it = FQtypes.begin(); it != FQtypes.end(); ++it) {    
688 >        int etid2 = Etids[ (*it) ];
689 >        int fqtid2 = FQtids[ (*it) ];
690 >        ElectrostaticAtomData eaData2 = ElectrostaticMap[ etid2 ];
691          RealType a = electrostaticAtomData.slaterZeta;
692          RealType b = eaData2.slaterZeta;
693          int m = electrostaticAtomData.slaterN;
694          int n = eaData2.slaterN;
695 <
695 >        
696          // Create the spline of the coulombic integral for s-type
697          // Slater orbitals.  Add a 2 angstrom safety window to deal
698          // with cutoffGroups that have charged atoms longer than the
699          // cutoffRadius away from each other.
700 <
700 >        
701          RealType rval;
702          RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
703          vector<RealType> rvals;
# Line 709 | Line 714 | namespace OpenMD {
714          
715          CubicSpline* J = new CubicSpline();
716          J->addPoints(rvals, Jvals);
717 <        
718 <        pair<AtomType*, AtomType*> key1, key2;
719 <        key1 = make_pair(atomType, atype2);
720 <        key2 = make_pair(atype2, atomType);
721 <        
717 <        Jij[key1] = J;
718 <        Jij[key2] = J;
719 <      }
720 <    }
721 <
717 >        Jij[fqtid][fqtid2] = J;
718 >        Jij[fqtid2].resize( nFlucq_ );
719 >        Jij[fqtid2][fqtid] = J;
720 >      }      
721 >    }      
722      return;
723    }
724    
# Line 744 | Line 744 | namespace OpenMD {
744  
745    void Electrostatic::calcForce(InteractionData &idat) {
746  
747 <    RealType C_a, C_b;  // Charges
748 <    Vector3d D_a, D_b;  // Dipoles (space-fixed)
749 <    Mat3x3d  Q_a, Q_b;  // Quadrupoles (space-fixed)
747 >    if (!initialized_) initialize();
748 >    
749 >    data1 = ElectrostaticMap[Etids[idat.atid1]];
750 >    data2 = ElectrostaticMap[Etids[idat.atid2]];
751  
752 <    RealType ri;                                 // Distance utility scalar
753 <    RealType rdDa, rdDb;                         // Dipole utility scalars
754 <    Vector3d rxDa, rxDb;                         // Dipole utility vectors
755 <    RealType rdQar, rdQbr, trQa, trQb;           // Quadrupole utility scalars
756 <    Vector3d Qar, Qbr, rQa, rQb, rxQar, rxQbr;   // Quadrupole utility vectors
757 <    RealType pref;
758 <
759 <    RealType DadDb, trQaQb, DadQbr, DbdQar;       // Cross-interaction scalars
759 <    RealType rQaQbr;
760 <    Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors
761 <    Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr;
762 <    Mat3x3d  QaQb;                                // Cross-interaction matrices
763 <
764 <    RealType U(0.0);  // Potential
765 <    Vector3d F(0.0);  // Force
766 <    Vector3d Ta(0.0); // Torque on site a
767 <    Vector3d Tb(0.0); // Torque on site b
768 <    Vector3d Ea(0.0); // Electric field at site a
769 <    Vector3d Eb(0.0); // Electric field at site b
770 <    RealType dUdCa(0.0); // fluctuating charge force at site a
771 <    RealType dUdCb(0.0); // fluctuating charge force at site a
752 >    U = 0.0;  // Potential
753 >    F.zero();  // Force
754 >    Ta.zero(); // Torque on site a
755 >    Tb.zero(); // Torque on site b
756 >    Ea.zero(); // Electric field at site a
757 >    Eb.zero(); // Electric field at site b
758 >    dUdCa = 0.0; // fluctuating charge force at site a
759 >    dUdCb = 0.0; // fluctuating charge force at site a
760      
761      // Indirect interactions mediated by the reaction field.
762 <    RealType indirect_Pot(0.0);  // Potential
763 <    Vector3d indirect_F(0.0);    // Force
764 <    Vector3d indirect_Ta(0.0);   // Torque on site a
765 <    Vector3d indirect_Tb(0.0);   // Torque on site b
762 >    indirect_Pot = 0.0;   // Potential
763 >    indirect_F.zero();    // Force
764 >    indirect_Ta.zero();   // Torque on site a
765 >    indirect_Tb.zero();   // Torque on site b
766  
767      // Excluded potential that is still computed for fluctuating charges
768 <    RealType excluded_Pot(0.0);
768 >    excluded_Pot= 0.0;
769  
782 –    RealType rfContrib, coulInt;
783 –    
784 –    // spline for coulomb integral
785 –    CubicSpline* J;
786 –
787 –    if (!initialized_) initialize();
788 –    
789 –    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
790 –    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
791 –    
770      // some variables we'll need independent of electrostatic type:
771  
772      ri = 1.0 /  *(idat.rij);
773 <    Vector3d rhat =  *(idat.d)  * ri;
773 >    rhat =  *(idat.d)  * ri;
774        
775      // logicals
776  
777 <    bool a_is_Charge = data1.is_Charge;
778 <    bool a_is_Dipole = data1.is_Dipole;
779 <    bool a_is_Quadrupole = data1.is_Quadrupole;
780 <    bool a_is_Fluctuating = data1.is_Fluctuating;
777 >    a_is_Charge = data1.is_Charge;
778 >    a_is_Dipole = data1.is_Dipole;
779 >    a_is_Quadrupole = data1.is_Quadrupole;
780 >    a_is_Fluctuating = data1.is_Fluctuating;
781  
782 <    bool b_is_Charge = data2.is_Charge;
783 <    bool b_is_Dipole = data2.is_Dipole;
784 <    bool b_is_Quadrupole = data2.is_Quadrupole;
785 <    bool b_is_Fluctuating = data2.is_Fluctuating;
782 >    b_is_Charge = data2.is_Charge;
783 >    b_is_Dipole = data2.is_Dipole;
784 >    b_is_Quadrupole = data2.is_Quadrupole;
785 >    b_is_Fluctuating = data2.is_Fluctuating;
786  
787      // Obtain all of the required radial function values from the
788      // spline structures:
# Line 909 | Line 887 | namespace OpenMD {
887          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
888                          + rdQbr * rhat * (dv22 - 2.0*v22or));
889      }
890 <    
890 >        
891 >
892      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
893 <      J = Jij[idat.atypes];
893 >      J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
894      }    
895 <    
895 >
896      if (a_is_Charge) {    
897        
898        if (b_is_Charge) {
899          pref =  pre11_ * *(idat.electroMult);      
900          U  += C_a * C_b * pref * v01;
901          F  += C_a * C_b * pref * dv01 * rhat;
902 <        
902 >
903          // If this is an excluded pair, there are still indirect
904          // interactions via the reaction field we must worry about:
905  
# Line 929 | Line 908 | namespace OpenMD {
908            indirect_Pot += rfContrib;
909            indirect_F   += rfContrib * 2.0 * ri * rhat;
910          }
911 <        
911 >
912          // Fluctuating charge forces are handled via Coulomb integrals
913          // for excluded pairs (i.e. those connected via bonds) and
914          // with the standard charge-charge interaction otherwise.
915  
916 <        if (idat.excluded) {          
916 >        if (idat.excluded) {
917            if (a_is_Fluctuating || b_is_Fluctuating) {
918              coulInt = J->getValueAt( *(idat.rij) );
919 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
920 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
921 <            excluded_Pot += C_a * C_b * coulInt;
943 <          }          
919 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
920 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
921 >          }
922          } else {
923            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
924 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
925 <        }
924 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
925 >        }              
926        }
927  
928        if (b_is_Dipole) {
# Line 1010 | Line 988 | namespace OpenMD {
988          F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
989          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
990          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
1013 –
991          // Even if we excluded this pair from direct interactions, we
992          // still have the reaction-field-mediated dipole-dipole
993          // interaction:
# Line 1070 | Line 1047 | namespace OpenMD {
1047          trQaQb = QaQb.trace();
1048          rQaQb = rhat * QaQb;
1049          QaQbr = QaQb * rhat;
1050 <        QaxQb = cross(Q_a, Q_b);
1050 >        QaxQb = mCross(Q_a, Q_b);
1051          rQaQbr = dot(rQa, Qbr);
1052          rQaxQbr = cross(rQa, Qbr);
1053          
# Line 1101 | Line 1078 | namespace OpenMD {
1078          //             + 4.0 * cross(rhat, QbQar)
1079  
1080          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1104 –
1105 –        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1081        }
1082      }
1083  
# Line 1149 | Line 1124 | namespace OpenMD {
1124  
1125      if (!initialized_) initialize();
1126  
1127 <    ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1127 >    ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]];
1128      
1129      // logicals
1130      bool i_is_Charge = data.is_Charge;
1131      bool i_is_Dipole = data.is_Dipole;
1132 +    bool i_is_Quadrupole = data.is_Quadrupole;
1133      bool i_is_Fluctuating = data.is_Fluctuating;
1134      RealType C_a = data.fixedCharge;  
1135 <    RealType self, preVal, DadDa;
1136 <    
1135 >    RealType self(0.0), preVal, DdD, trQ, trQQ;
1136 >
1137 >    if (i_is_Dipole) {
1138 >      DdD = data.dipole.lengthSquare();
1139 >    }
1140 >        
1141      if (i_is_Fluctuating) {
1142        C_a += *(sdat.flucQ);
1143 <      // dVdFQ is really a force, so this is negative the derivative
1164 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1143 >      *(sdat.flucQfrc) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1144        (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1145          (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1146      }
# Line 1178 | Line 1157 | namespace OpenMD {
1157        }
1158  
1159        if (i_is_Dipole) {
1160 <        DadDa = data.dipole.lengthSquare();
1182 <        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DadDa;
1160 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DdD;
1161        }
1162        
1163        break;
1164        
1165      case esm_SHIFTED_FORCE:
1166      case esm_SHIFTED_POTENTIAL:
1167 <      if (i_is_Charge) {
1168 <        self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1169 <        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1167 >    case esm_TAYLOR_SHIFTED:
1168 >    case esm_EWALD_FULL:
1169 >      if (i_is_Charge)
1170 >        self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge));      
1171 >      if (i_is_Dipole)
1172 >        self += selfMult2_ * pre22_ * DdD;      
1173 >      if (i_is_Quadrupole) {
1174 >        trQ = data.quadrupole.trace();
1175 >        trQQ = (data.quadrupole * data.quadrupole).trace();
1176 >        self += selfMult4_ * pre44_ * (2.0*trQQ + trQ*trQ);
1177 >        if (i_is_Charge)
1178 >          self -= selfMult2_ * pre14_ * 2.0 * C_a * trQ;
1179        }
1180 +      (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;      
1181        break;
1182      default:
1183        break;
# Line 1203 | Line 1191 | namespace OpenMD {
1191      // cases.
1192      return 12.0;
1193    }
1194 +
1195 +
1196 +  void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1197 +    
1198 +    RealType kPot = 0.0;
1199 +    RealType kVir = 0.0;
1200 +    
1201 +    const RealType mPoleConverter = 0.20819434; // converts from the
1202 +                                                // internal units of
1203 +                                                // Debye (for dipoles)
1204 +                                                // or Debye-angstroms
1205 +                                                // (for quadrupoles) to
1206 +                                                // electron angstroms or
1207 +                                                // electron-angstroms^2
1208 +    
1209 +    const RealType eConverter = 332.0637778; // convert the
1210 +                                             // Charge-Charge
1211 +                                             // electrostatic
1212 +                                             // interactions into kcal /
1213 +                                             // mol assuming distances
1214 +                                             // are measured in
1215 +                                             // angstroms.
1216 +
1217 +    Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat();
1218 +    Vector3d box = hmat.diagonals();
1219 +    RealType boxMax = box.max();
1220 +    
1221 +    //int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) );
1222 +    int kMax = 7;
1223 +    int kSqMax = kMax*kMax + 2;
1224 +    
1225 +    int kLimit = kMax+1;
1226 +    int kLim2 = 2*kMax+1;
1227 +    int kSqLim = kSqMax;
1228 +    
1229 +    vector<RealType> AK(kSqLim+1, 0.0);
1230 +    RealType xcl = 2.0 * M_PI / box.x();
1231 +    RealType ycl = 2.0 * M_PI / box.y();
1232 +    RealType zcl = 2.0 * M_PI / box.z();
1233 +    RealType rcl = 2.0 * M_PI / boxMax;
1234 +    RealType rvol = 2.0 * M_PI /(box.x() * box.y() * box.z());
1235 +    
1236 +    if(dampingAlpha_ < 1.0e-12) return;
1237 +    
1238 +    RealType ralph = -0.25/pow(dampingAlpha_,2);
1239 +    
1240 +    // Calculate and store exponential factors  
1241 +    
1242 +    vector<vector<RealType> > elc;
1243 +    vector<vector<RealType> > emc;
1244 +    vector<vector<RealType> > enc;
1245 +    vector<vector<RealType> > els;
1246 +    vector<vector<RealType> > ems;
1247 +    vector<vector<RealType> > ens;
1248 +
1249 +    
1250 +    int nMax = info_->getNAtoms();
1251 +    
1252 +    elc.resize(kLimit+1);
1253 +    emc.resize(kLimit+1);
1254 +    enc.resize(kLimit+1);
1255 +    els.resize(kLimit+1);
1256 +    ems.resize(kLimit+1);
1257 +    ens.resize(kLimit+1);
1258 +
1259 +    for (int j = 0; j < kLimit+1; j++) {
1260 +      elc[j].resize(nMax);
1261 +      emc[j].resize(nMax);
1262 +      enc[j].resize(nMax);
1263 +      els[j].resize(nMax);
1264 +      ems[j].resize(nMax);
1265 +      ens[j].resize(nMax);
1266 +    }
1267 +    
1268 +    Vector3d t( 2.0 * M_PI );
1269 +    t.Vdiv(t, box);
1270 +
1271 +    
1272 +    SimInfo::MoleculeIterator mi;
1273 +    Molecule::AtomIterator ai;
1274 +    int i;
1275 +    Vector3d r;
1276 +    Vector3d tt;
1277 +    
1278 +    for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1279 +         mol = info_->nextMolecule(mi)) {
1280 +      for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1281 +          atom = mol->nextAtom(ai)) {  
1282 +        
1283 +        i = atom->getLocalIndex();
1284 +        r = atom->getPos();
1285 +        info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
1286 +        
1287 +        tt.Vmul(t, r);
1288 +
1289 +        elc[1][i] = 1.0;
1290 +        emc[1][i] = 1.0;
1291 +        enc[1][i] = 1.0;
1292 +        els[1][i] = 0.0;
1293 +        ems[1][i] = 0.0;
1294 +        ens[1][i] = 0.0;
1295 +
1296 +        elc[2][i] = cos(tt.x());
1297 +        emc[2][i] = cos(tt.y());
1298 +        enc[2][i] = cos(tt.z());
1299 +        els[2][i] = sin(tt.x());
1300 +        ems[2][i] = sin(tt.y());
1301 +        ens[2][i] = sin(tt.z());
1302 +        
1303 +        for(int l = 3; l <= kLimit; l++) {
1304 +          elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1305 +          emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1306 +          enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1307 +          els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1308 +          ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1309 +          ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1310 +        }
1311 +      }
1312 +    }
1313 +    
1314 +    // Calculate and store AK coefficients:
1315 +    
1316 +    RealType eksq = 1.0;
1317 +    RealType expf = 0.0;
1318 +    if (ralph < 0.0) expf = exp(ralph*rcl*rcl);
1319 +    for (i = 1; i <= kSqLim; i++) {
1320 +      RealType rksq = float(i)*rcl*rcl;
1321 +      eksq = expf*eksq;
1322 +      AK[i] = eConverter * eksq/rksq;
1323 +    }
1324 +    
1325 +    /*
1326 +     * Loop over all k vectors k = 2 pi (ll/Lx, mm/Ly, nn/Lz)
1327 +     * the values of ll, mm and nn are selected so that the symmetry of
1328 +     * reciprocal lattice is taken into account i.e. the following
1329 +     * rules apply.
1330 +     *
1331 +     * ll ranges over the values 0 to kMax only.
1332 +     *
1333 +     * mm ranges over 0 to kMax when ll=0 and over
1334 +     *            -kMax to kMax otherwise.
1335 +     * nn ranges over 1 to kMax when ll=mm=0 and over
1336 +     *            -kMax to kMax otherwise.
1337 +     *
1338 +     * Hence the result of the summation must be doubled at the end.    
1339 +     */
1340 +    
1341 +    std::vector<RealType> clm(nMax, 0.0);
1342 +    std::vector<RealType> slm(nMax, 0.0);
1343 +    std::vector<RealType> ckr(nMax, 0.0);
1344 +    std::vector<RealType> skr(nMax, 0.0);
1345 +    std::vector<RealType> ckc(nMax, 0.0);
1346 +    std::vector<RealType> cks(nMax, 0.0);
1347 +    std::vector<RealType> dkc(nMax, 0.0);
1348 +    std::vector<RealType> dks(nMax, 0.0);
1349 +    std::vector<RealType> qkc(nMax, 0.0);
1350 +    std::vector<RealType> qks(nMax, 0.0);
1351 +    std::vector<Vector3d> dxk(nMax, V3Zero);
1352 +    std::vector<Vector3d> qxk(nMax, V3Zero);
1353 +    RealType rl, rm, rn;
1354 +    Vector3d kVec;
1355 +    Vector3d Qk;
1356 +    Mat3x3d k2;
1357 +    RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1358 +    int atid;
1359 +    ElectrostaticAtomData data;
1360 +    RealType C, dk, qk;
1361 +    Vector3d D;
1362 +    Mat3x3d  Q;
1363 +
1364 +    int mMin = kLimit;
1365 +    int nMin = kLimit + 1;
1366 +    for (int l = 1; l <= kLimit; l++) {
1367 +      int ll = l - 1;
1368 +      rl = xcl * float(ll);
1369 +      for (int mmm = mMin; mmm <= kLim2; mmm++) {
1370 +        int mm = mmm - kLimit;
1371 +        int m = abs(mm) + 1;
1372 +        rm = ycl * float(mm);
1373 +        // Set temporary products of exponential terms
1374 +        for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1375 +             mol = info_->nextMolecule(mi)) {
1376 +          for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1377 +              atom = mol->nextAtom(ai)) {
1378 +            
1379 +            i = atom->getLocalIndex();
1380 +            if(mm < 0) {
1381 +              clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1382 +              slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1383 +            } else {
1384 +              clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1385 +              slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1386 +            }
1387 +          }
1388 +        }
1389 +        for (int nnn = nMin; nnn <= kLim2; nnn++) {
1390 +          int nn = nnn - kLimit;          
1391 +          int n = abs(nn) + 1;
1392 +          rn = zcl * float(nn);
1393 +          // Test on magnitude of k vector:
1394 +          int kk=ll*ll + mm*mm + nn*nn;
1395 +          if(kk <= kSqLim) {
1396 +            kVec = Vector3d(rl, rm, rn);
1397 +            k2 = outProduct(kVec, kVec);
1398 +            // Calculate exp(ikr) terms
1399 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1400 +                 mol = info_->nextMolecule(mi)) {
1401 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1402 +                  atom = mol->nextAtom(ai)) {
1403 +                i = atom->getLocalIndex();
1404 +                
1405 +                if (nn < 0) {
1406 +                  ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1407 +                  skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1408 +
1409 +                } else {
1410 +                  ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1411 +                  skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1412 +                }
1413 +              }
1414 +            }
1415 +            
1416 +            // Calculate scalar and vector products for each site:
1417 +            
1418 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1419 +                 mol = info_->nextMolecule(mi)) {
1420 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1421 +                  atom = mol->nextAtom(ai)) {
1422 +                i = atom->getLocalIndex();
1423 +                int atid = atom->getAtomType()->getIdent();
1424 +                data = ElectrostaticMap[Etids[atid]];
1425 +                              
1426 +                if (data.is_Charge) {
1427 +                  C = data.fixedCharge;
1428 +                  if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1429 +                  ckc[i] = C * ckr[i];
1430 +                  cks[i] = C * skr[i];
1431 +                }
1432 +                
1433 +                if (data.is_Dipole) {
1434 +                  D = atom->getDipole() * mPoleConverter;
1435 +                  dk = dot(D, kVec);
1436 +                  dxk[i] = cross(D, kVec);
1437 +                  dkc[i] = dk * ckr[i];
1438 +                  dks[i] = dk * skr[i];
1439 +                }
1440 +                if (data.is_Quadrupole) {
1441 +                  Q = atom->getQuadrupole() * mPoleConverter;
1442 +                  Qk = Q * kVec;                  
1443 +                  qk = dot(kVec, Qk);
1444 +                  qxk[i] = -cross(kVec, Qk);
1445 +                  qkc[i] = qk * ckr[i];
1446 +                  qks[i] = qk * skr[i];
1447 +                }              
1448 +              }
1449 +            }
1450 +
1451 +            // calculate vector sums
1452 +            
1453 +            ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1454 +            ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1455 +            dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1456 +            dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1457 +            qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1458 +            qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1459 +            
1460 + #ifdef IS_MPI
1461 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
1462 +                                      MPI::SUM);
1463 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE,
1464 +                                      MPI::SUM);
1465 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE,
1466 +                                      MPI::SUM);
1467 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE,
1468 +                                      MPI::SUM);
1469 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE,
1470 +                                      MPI::SUM);
1471 +            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE,
1472 +                                      MPI::SUM);
1473 + #endif        
1474 +            
1475 +            // Accumulate potential energy and virial contribution:
1476 +
1477 +            kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1478 +                                         + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1479 +
1480 +            kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1481 +                                          +4.0*(ckss*dkcs-ckcs*dkss)
1482 +                                          +3.0*(dkcs*dkcs+dkss*dkss)
1483 +                                          -6.0*(ckss*qkss+ckcs*qkcs)
1484 +                                          +8.0*(dkss*qkcs-dkcs*qkss)
1485 +                                          +5.0*(qkss*qkss+qkcs*qkcs));
1486 +            
1487 +            // Calculate force and torque for each site:
1488 +            
1489 +            for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1490 +                 mol = info_->nextMolecule(mi)) {
1491 +              for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1492 +                  atom = mol->nextAtom(ai)) {
1493 +                
1494 +                i = atom->getLocalIndex();
1495 +                atid = atom->getAtomType()->getIdent();
1496 +                data = ElectrostaticMap[Etids[atid]];
1497 +
1498 +                RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1499 +                                     - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1500 +                RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1501 +                                         -ckr[i]*(ckss+dkcs-qkss));
1502 +                RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1503 +                                            +skr[i]*(ckss+dkcs-qkss));
1504 +              
1505 +                atom->addFrc( 4.0 * rvol * qfrc * kVec );
1506 +                
1507 +                if (data.is_Dipole) {
1508 +                  atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1509 +                }
1510 +                if (data.is_Quadrupole) {
1511 +                  atom->addTrq( 4.0 * rvol * qtrq2 * qxk[i] );
1512 +                }
1513 +              }
1514 +            }
1515 +          }
1516 +        }
1517 +        nMin = 1;
1518 +      }
1519 +      mMin = 1;
1520 +    }
1521 +    pot += kPot;  
1522 +  }
1523   }

Diff Legend

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