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 1895 by gezelter, Mon Jul 1 21:09:37 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 398 | Line 410 | namespace OpenMD {
410          v11 = g - gc - rmRc*hc;
411          v21 = g*ri - gc*ric - rmRc*(hc - gc*ric)*ric;
412          v22 = h - g*ri - (hc - gc*ric) - rmRc*(sc - (hc - gc*ric)*ric);
413 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
413 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
414          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric)
415            - rmRc*(tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
416          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2
# Line 455 | Line 467 | namespace OpenMD {
467          v11 = g - gc;
468          v21 = g*ri - gc*ric;
469          v22 = h - g*ri - (hc - gc*ric);
470 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
470 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric;
471          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
472          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
473          v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;        
# Line 515 | Line 527 | namespace OpenMD {
527          v11 = g - gc;
528          v21 = g*ri - gc*ric;
529          v22 = h - g*ri - (hc - gc*ric);
530 <        v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
530 >        v31 = (h-g*ri)*ri - (hc-gc*ric)*ric;
531          v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
532          v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
533          v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;        
# Line 625 | Line 637 | namespace OpenMD {
637    }
638        
639    void Electrostatic::addType(AtomType* atomType){
640 <
640 >    
641      ElectrostaticAtomData electrostaticAtomData;
642      electrostaticAtomData.is_Charge = false;
643      electrostaticAtomData.is_Dipole = false;
# Line 661 | Line 673 | namespace OpenMD {
673        electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
674      }
675  
676 <    pair<map<int,AtomType*>::iterator,bool> ret;    
677 <    ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
678 <                                                        atomType) );
676 >    int atid = atomType->getIdent();
677 >    int etid = Etypes.size();
678 >    int fqtid = FQtypes.size();
679 >
680 >    pair<set<int>::iterator,bool> ret;    
681 >    ret = Etypes.insert( atid );
682      if (ret.second == false) {
683        sprintf( painCave.errMsg,
684                 "Electrostatic already had a previous entry with ident %d\n",
685 <               atomType->getIdent() );
685 >               atid);
686        painCave.severity = OPENMD_INFO;
687        painCave.isFatal = 0;
688        simError();        
689      }
690      
691 <    ElectrostaticMap[atomType] = electrostaticAtomData;  
691 >    Etids[ atid ] = etid;
692 >    ElectrostaticMap.push_back(electrostaticAtomData);
693  
694 <    // Now, iterate over all known types and add to the mixing map:
695 <    
696 <    map<AtomType*, ElectrostaticAtomData>::iterator it;
697 <    for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
698 <      AtomType* atype2 = (*it).first;
699 <      ElectrostaticAtomData eaData2 = (*it).second;
700 <      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
701 <        
694 >    if (electrostaticAtomData.is_Fluctuating) {
695 >      ret = FQtypes.insert( atid );
696 >      if (ret.second == false) {
697 >        sprintf( painCave.errMsg,
698 >                 "Electrostatic already had a previous fluctuating charge entry with ident %d\n",
699 >                 atid );
700 >        painCave.severity = OPENMD_INFO;
701 >        painCave.isFatal = 0;
702 >        simError();        
703 >      }
704 >      FQtids[atid] = fqtid;
705 >      Jij[fqtid].resize(nFlucq_);
706 >
707 >      // Now, iterate over all known fluctuating and add to the coulomb integral map:
708 >      
709 >      std::set<int>::iterator it;
710 >      for( it = FQtypes.begin(); it != FQtypes.end(); ++it) {    
711 >        int etid2 = Etids[ (*it) ];
712 >        int fqtid2 = FQtids[ (*it) ];
713 >        ElectrostaticAtomData eaData2 = ElectrostaticMap[ etid2 ];
714          RealType a = electrostaticAtomData.slaterZeta;
715          RealType b = eaData2.slaterZeta;
716          int m = electrostaticAtomData.slaterN;
717          int n = eaData2.slaterN;
718 <
718 >        
719          // Create the spline of the coulombic integral for s-type
720          // Slater orbitals.  Add a 2 angstrom safety window to deal
721          // with cutoffGroups that have charged atoms longer than the
722          // cutoffRadius away from each other.
723 <
723 >        
724          RealType rval;
725          RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
726          vector<RealType> rvals;
# Line 709 | Line 737 | namespace OpenMD {
737          
738          CubicSpline* J = new CubicSpline();
739          J->addPoints(rvals, Jvals);
740 <        
741 <        pair<AtomType*, AtomType*> key1, key2;
742 <        key1 = make_pair(atomType, atype2);
743 <        key2 = make_pair(atype2, atomType);
744 <        
717 <        Jij[key1] = J;
718 <        Jij[key2] = J;
719 <      }
720 <    }
721 <
740 >        Jij[fqtid][fqtid2] = J;
741 >        Jij[fqtid2].resize( nFlucq_ );
742 >        Jij[fqtid2][fqtid] = J;
743 >      }      
744 >    }      
745      return;
746    }
747    
# Line 744 | Line 767 | namespace OpenMD {
767  
768    void Electrostatic::calcForce(InteractionData &idat) {
769  
770 <    RealType C_a, C_b;  // Charges
771 <    Vector3d D_a, D_b;  // Dipoles (space-fixed)
772 <    Mat3x3d  Q_a, Q_b;  // Quadrupoles (space-fixed)
770 >    if (!initialized_) initialize();
771 >    
772 >    data1 = ElectrostaticMap[Etids[idat.atid1]];
773 >    data2 = ElectrostaticMap[Etids[idat.atid2]];
774  
775 <    RealType ri;                                 // Distance utility scalar
776 <    RealType rdDa, rdDb;                         // Dipole utility scalars
777 <    Vector3d rxDa, rxDb;                         // Dipole utility vectors
778 <    RealType rdQar, rdQbr, trQa, trQb;           // Quadrupole utility scalars
779 <    Vector3d Qar, Qbr, rQa, rQb, rxQar, rxQbr;   // Quadrupole utility vectors
780 <    RealType pref;
781 <
782 <    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
775 >    U = 0.0;  // Potential
776 >    F.zero();  // Force
777 >    Ta.zero(); // Torque on site a
778 >    Tb.zero(); // Torque on site b
779 >    Ea.zero(); // Electric field at site a
780 >    Eb.zero(); // Electric field at site b
781 >    dUdCa = 0.0; // fluctuating charge force at site a
782 >    dUdCb = 0.0; // fluctuating charge force at site a
783      
784      // Indirect interactions mediated by the reaction field.
785 <    RealType indirect_Pot(0.0);  // Potential
786 <    Vector3d indirect_F(0.0);    // Force
787 <    Vector3d indirect_Ta(0.0);   // Torque on site a
788 <    Vector3d indirect_Tb(0.0);   // Torque on site b
785 >    indirect_Pot = 0.0;   // Potential
786 >    indirect_F.zero();    // Force
787 >    indirect_Ta.zero();   // Torque on site a
788 >    indirect_Tb.zero();   // Torque on site b
789  
790      // Excluded potential that is still computed for fluctuating charges
791 <    RealType excluded_Pot(0.0);
791 >    excluded_Pot= 0.0;
792  
782 –    RealType rfContrib, coulInt;
783 –    
784 –    // spline for coulomb integral
785 –    CubicSpline* J;
793  
787 –    if (!initialized_) initialize();
788 –    
789 –    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
790 –    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
791 –    
794      // some variables we'll need independent of electrostatic type:
795  
796      ri = 1.0 /  *(idat.rij);
797 <    Vector3d rhat =  *(idat.d)  * ri;
797 >    rhat =  *(idat.d)  * ri;
798        
799      // logicals
800  
801 <    bool a_is_Charge = data1.is_Charge;
802 <    bool a_is_Dipole = data1.is_Dipole;
803 <    bool a_is_Quadrupole = data1.is_Quadrupole;
804 <    bool a_is_Fluctuating = data1.is_Fluctuating;
801 >    a_is_Charge = data1.is_Charge;
802 >    a_is_Dipole = data1.is_Dipole;
803 >    a_is_Quadrupole = data1.is_Quadrupole;
804 >    a_is_Fluctuating = data1.is_Fluctuating;
805  
806 <    bool b_is_Charge = data2.is_Charge;
807 <    bool b_is_Dipole = data2.is_Dipole;
808 <    bool b_is_Quadrupole = data2.is_Quadrupole;
809 <    bool b_is_Fluctuating = data2.is_Fluctuating;
806 >    b_is_Charge = data2.is_Charge;
807 >    b_is_Dipole = data2.is_Dipole;
808 >    b_is_Quadrupole = data2.is_Quadrupole;
809 >    b_is_Fluctuating = data2.is_Fluctuating;
810  
811      // Obtain all of the required radial function values from the
812      // spline structures:
# Line 911 | Line 913 | namespace OpenMD {
913      }
914      
915      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
916 <      J = Jij[idat.atypes];
916 >      J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
917      }    
918      
919      if (a_is_Charge) {    
# Line 1149 | Line 1151 | namespace OpenMD {
1151  
1152      if (!initialized_) initialize();
1153  
1154 <    ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1154 >    ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]];
1155      
1156      // logicals
1157      bool i_is_Charge = data.is_Charge;

Diff Legend

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