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 1907 by gezelter, Fri Jul 19 18:18:27 2013 UTC

# Line 213 | Line 213 | namespace OpenMD {
213        haveDampingAlpha_ = true;
214      }
215  
216 <    // find all of the Electrostatic atom Types:
217 <    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
218 <    ForceField::AtomTypeContainer::MapTypeIterator i;
219 <    AtomType* at;
216 >    Etypes.clear();
217 >    Etids.clear();
218 >    FQtypes.clear();
219 >    FQtids.clear();
220 >    ElectrostaticMap.clear();
221 >    Jij.clear();
222 >    nElectro_ = 0;
223 >    nFlucq_ = 0;
224 >
225 >    Etids.resize( forceField_->getNAtomType(), -1);
226 >    FQtids.resize( forceField_->getNAtomType(), -1);
227 >
228 >    set<AtomType*>::iterator at;
229 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {    
230 >      if ((*at)->isElectrostatic()) nElectro_++;
231 >      if ((*at)->isFluctuatingCharge()) nFlucq_++;
232 >    }
233      
234 <    for (at = atomTypes->beginType(i); at != NULL;
235 <         at = atomTypes->nextType(i)) {
236 <      
237 <      if (at->isElectrostatic())
225 <        addType(at);
234 >    Jij.resize(nFlucq_);
235 >
236 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
237 >      if ((*at)->isElectrostatic()) addType(*at);
238      }  
239      
240      if (summationMethod_ == esm_REACTION_FIELD) {
# Line 250 | Line 262 | namespace OpenMD {
262        b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
263        b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
264        b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
265 <      selfMult_ = b0c + a2 * invArootPi;
265 >      //selfMult1_ = - 2.0 * a2 * invArootPi;
266 >      //selfMult2_ = - 4.0 * a2 * a2 * invArootPi / 3.0;
267 >      //selfMult4_ = - 8.0 * a2 * a2 * a2 * invArootPi / 5.0;
268 >      // Half the Smith self piece:
269 >      selfMult1_ = - a2 * invArootPi;
270 >      selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0;
271 >      selfMult4_ = - 4.0 * a2 * a2 * a2 * invArootPi / 5.0;
272      } else {
273        a2 = 0.0;
274        b0c = 1.0 / r;
# Line 259 | Line 277 | namespace OpenMD {
277        b3c = (5.0 * b2c) / r2;
278        b4c = (7.0 * b3c) / r2;
279        b5c = (9.0 * b4c) / r2;
280 <      selfMult_ = b0c;
280 >      selfMult1_ = 0.0;
281 >      selfMult2_ = 0.0;
282 >      selfMult4_ = 0.0;
283      }
284  
285      // higher derivatives of B_0 at the cutoff radius:
# Line 267 | Line 287 | namespace OpenMD {
287      db0c_2 =     -b1c + r2 * b2c;
288      db0c_3 =          3.0*r*b2c  - r2*r*b3c;
289      db0c_4 =          3.0*b2c  - 6.0*r2*b3c     + r2*r2*b4c;
290 <    db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;
290 >    db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;  
291      
292 +    selfMult1_ -= b0c;
293 +    selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) /  3.0;
294 +    selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0;
295  
296      // working variables for the splines:
297      RealType ri, ri2;
# Line 398 | Line 421 | namespace OpenMD {
421          v11 = g - gc - rmRc*hc;
422          v21 = g*ri - gc*ric - rmRc*(hc - gc*ric)*ric;
423          v22 = h - g*ri - (hc - gc*ric) - rmRc*(sc - (hc - gc*ric)*ric);
424 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
424 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
425          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric)
426            - rmRc*(tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
427          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2
# Line 455 | Line 478 | namespace OpenMD {
478          v11 = g - gc;
479          v21 = g*ri - gc*ric;
480          v22 = h - g*ri - (hc - gc*ric);
481 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
481 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric;
482          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
483          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
484          v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;        
# Line 515 | Line 538 | namespace OpenMD {
538          v11 = g - gc;
539          v21 = g*ri - gc*ric;
540          v22 = h - g*ri - (hc - gc*ric);
541 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
541 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric;
542          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
543          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
544          v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;        
# Line 625 | Line 648 | namespace OpenMD {
648    }
649        
650    void Electrostatic::addType(AtomType* atomType){
651 <
651 >    
652      ElectrostaticAtomData electrostaticAtomData;
653      electrostaticAtomData.is_Charge = false;
654      electrostaticAtomData.is_Dipole = false;
# Line 661 | Line 684 | namespace OpenMD {
684        electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
685      }
686  
687 <    pair<map<int,AtomType*>::iterator,bool> ret;    
688 <    ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
689 <                                                        atomType) );
687 >    int atid = atomType->getIdent();
688 >    int etid = Etypes.size();
689 >    int fqtid = FQtypes.size();
690 >
691 >    pair<set<int>::iterator,bool> ret;    
692 >    ret = Etypes.insert( atid );
693      if (ret.second == false) {
694        sprintf( painCave.errMsg,
695                 "Electrostatic already had a previous entry with ident %d\n",
696 <               atomType->getIdent() );
696 >               atid);
697        painCave.severity = OPENMD_INFO;
698        painCave.isFatal = 0;
699        simError();        
700      }
701      
702 <    ElectrostaticMap[atomType] = electrostaticAtomData;  
702 >    Etids[ atid ] = etid;
703 >    ElectrostaticMap.push_back(electrostaticAtomData);
704  
705 <    // Now, iterate over all known types and add to the mixing map:
706 <    
707 <    map<AtomType*, ElectrostaticAtomData>::iterator it;
708 <    for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
709 <      AtomType* atype2 = (*it).first;
710 <      ElectrostaticAtomData eaData2 = (*it).second;
711 <      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
712 <        
705 >    if (electrostaticAtomData.is_Fluctuating) {
706 >      ret = FQtypes.insert( atid );
707 >      if (ret.second == false) {
708 >        sprintf( painCave.errMsg,
709 >                 "Electrostatic already had a previous fluctuating charge entry with ident %d\n",
710 >                 atid );
711 >        painCave.severity = OPENMD_INFO;
712 >        painCave.isFatal = 0;
713 >        simError();        
714 >      }
715 >      FQtids[atid] = fqtid;
716 >      Jij[fqtid].resize(nFlucq_);
717 >
718 >      // Now, iterate over all known fluctuating and add to the coulomb integral map:
719 >      
720 >      std::set<int>::iterator it;
721 >      for( it = FQtypes.begin(); it != FQtypes.end(); ++it) {    
722 >        int etid2 = Etids[ (*it) ];
723 >        int fqtid2 = FQtids[ (*it) ];
724 >        ElectrostaticAtomData eaData2 = ElectrostaticMap[ etid2 ];
725          RealType a = electrostaticAtomData.slaterZeta;
726          RealType b = eaData2.slaterZeta;
727          int m = electrostaticAtomData.slaterN;
728          int n = eaData2.slaterN;
729 <
729 >        
730          // Create the spline of the coulombic integral for s-type
731          // Slater orbitals.  Add a 2 angstrom safety window to deal
732          // with cutoffGroups that have charged atoms longer than the
733          // cutoffRadius away from each other.
734 <
734 >        
735          RealType rval;
736          RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
737          vector<RealType> rvals;
# Line 709 | Line 748 | namespace OpenMD {
748          
749          CubicSpline* J = new CubicSpline();
750          J->addPoints(rvals, Jvals);
751 <        
752 <        pair<AtomType*, AtomType*> key1, key2;
753 <        key1 = make_pair(atomType, atype2);
754 <        key2 = make_pair(atype2, atomType);
755 <        
717 <        Jij[key1] = J;
718 <        Jij[key2] = J;
719 <      }
720 <    }
721 <
751 >        Jij[fqtid][fqtid2] = J;
752 >        Jij[fqtid2].resize( nFlucq_ );
753 >        Jij[fqtid2][fqtid] = J;
754 >      }      
755 >    }      
756      return;
757    }
758    
# Line 744 | Line 778 | namespace OpenMD {
778  
779    void Electrostatic::calcForce(InteractionData &idat) {
780  
781 <    RealType C_a, C_b;  // Charges
782 <    Vector3d D_a, D_b;  // Dipoles (space-fixed)
783 <    Mat3x3d  Q_a, Q_b;  // Quadrupoles (space-fixed)
781 >    if (!initialized_) initialize();
782 >    
783 >    data1 = ElectrostaticMap[Etids[idat.atid1]];
784 >    data2 = ElectrostaticMap[Etids[idat.atid2]];
785  
786 <    RealType ri;                                 // Distance utility scalar
787 <    RealType rdDa, rdDb;                         // Dipole utility scalars
788 <    Vector3d rxDa, rxDb;                         // Dipole utility vectors
789 <    RealType rdQar, rdQbr, trQa, trQb;           // Quadrupole utility scalars
790 <    Vector3d Qar, Qbr, rQa, rQb, rxQar, rxQbr;   // Quadrupole utility vectors
791 <    RealType pref;
792 <
793 <    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
786 >    U = 0.0;  // Potential
787 >    F.zero();  // Force
788 >    Ta.zero(); // Torque on site a
789 >    Tb.zero(); // Torque on site b
790 >    Ea.zero(); // Electric field at site a
791 >    Eb.zero(); // Electric field at site b
792 >    dUdCa = 0.0; // fluctuating charge force at site a
793 >    dUdCb = 0.0; // fluctuating charge force at site a
794      
795      // Indirect interactions mediated by the reaction field.
796 <    RealType indirect_Pot(0.0);  // Potential
797 <    Vector3d indirect_F(0.0);    // Force
798 <    Vector3d indirect_Ta(0.0);   // Torque on site a
799 <    Vector3d indirect_Tb(0.0);   // Torque on site b
796 >    indirect_Pot = 0.0;   // Potential
797 >    indirect_F.zero();    // Force
798 >    indirect_Ta.zero();   // Torque on site a
799 >    indirect_Tb.zero();   // Torque on site b
800  
801      // Excluded potential that is still computed for fluctuating charges
802 <    RealType excluded_Pot(0.0);
802 >    excluded_Pot= 0.0;
803  
782 –    RealType rfContrib, coulInt;
783 –    
784 –    // spline for coulomb integral
785 –    CubicSpline* J;
804  
787 –    if (!initialized_) initialize();
788 –    
789 –    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
790 –    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
791 –    
805      // some variables we'll need independent of electrostatic type:
806  
807      ri = 1.0 /  *(idat.rij);
808 <    Vector3d rhat =  *(idat.d)  * ri;
808 >    rhat =  *(idat.d)  * ri;
809        
810      // logicals
811  
812 <    bool a_is_Charge = data1.is_Charge;
813 <    bool a_is_Dipole = data1.is_Dipole;
814 <    bool a_is_Quadrupole = data1.is_Quadrupole;
815 <    bool a_is_Fluctuating = data1.is_Fluctuating;
812 >    a_is_Charge = data1.is_Charge;
813 >    a_is_Dipole = data1.is_Dipole;
814 >    a_is_Quadrupole = data1.is_Quadrupole;
815 >    a_is_Fluctuating = data1.is_Fluctuating;
816  
817 <    bool b_is_Charge = data2.is_Charge;
818 <    bool b_is_Dipole = data2.is_Dipole;
819 <    bool b_is_Quadrupole = data2.is_Quadrupole;
820 <    bool b_is_Fluctuating = data2.is_Fluctuating;
817 >    b_is_Charge = data2.is_Charge;
818 >    b_is_Dipole = data2.is_Dipole;
819 >    b_is_Quadrupole = data2.is_Quadrupole;
820 >    b_is_Fluctuating = data2.is_Fluctuating;
821  
822      // Obtain all of the required radial function values from the
823      // spline structures:
# Line 911 | Line 924 | namespace OpenMD {
924      }
925      
926      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
927 <      J = Jij[idat.atypes];
927 >      J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
928      }    
929      
930      if (a_is_Charge) {    
# Line 1102 | Line 1115 | namespace OpenMD {
1115  
1116          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1117  
1105 –        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1118        }
1119      }
1120  
# Line 1149 | Line 1161 | namespace OpenMD {
1161  
1162      if (!initialized_) initialize();
1163  
1164 <    ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1164 >    ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]];
1165      
1166      // logicals
1167      bool i_is_Charge = data.is_Charge;
1168      bool i_is_Dipole = data.is_Dipole;
1169 +    bool i_is_Quadrupole = data.is_Quadrupole;
1170      bool i_is_Fluctuating = data.is_Fluctuating;
1171      RealType C_a = data.fixedCharge;  
1172 <    RealType self, preVal, DadDa;
1173 <    
1172 >    RealType self(0.0), preVal, DdD, trQ, trQQ;
1173 >
1174 >    if (i_is_Dipole) {
1175 >      DdD = data.dipole.lengthSquare();
1176 >    }
1177 >        
1178      if (i_is_Fluctuating) {
1179        C_a += *(sdat.flucQ);
1180        // dVdFQ is really a force, so this is negative the derivative
# Line 1178 | Line 1195 | namespace OpenMD {
1195        }
1196  
1197        if (i_is_Dipole) {
1198 <        DadDa = data.dipole.lengthSquare();
1182 <        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DadDa;
1198 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DdD;
1199        }
1200        
1201        break;
1202        
1203      case esm_SHIFTED_FORCE:
1204      case esm_SHIFTED_POTENTIAL:
1205 <      if (i_is_Charge) {
1206 <        self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1207 <        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1205 >    case esm_TAYLOR_SHIFTED:
1206 >      if (i_is_Charge)
1207 >        self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge));      
1208 >      if (i_is_Dipole)
1209 >        self += selfMult2_ * pre22_ * DdD;      
1210 >      if (i_is_Quadrupole) {
1211 >        trQ = data.quadrupole.trace();
1212 >        trQQ = (data.quadrupole * data.quadrupole).trace();
1213 >        self += selfMult4_ * pre44_ * (2.0*trQQ + trQ*trQ);
1214 >        if (i_is_Charge)
1215 >          self -= selfMult2_ * pre14_ * 2.0 * C_a * trQ;
1216        }
1217 +      (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;      
1218        break;
1219      default:
1220        break;

Diff Legend

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