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 1925 by gezelter, Wed Aug 7 15:24:16 2013 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 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  
# Line 57 | Line 61
61   #include "math/erfc.hpp"
62   #include "math/SquareMatrix.hpp"
63   #include "primitives/Molecule.hpp"
64 < #ifdef IS_MPI
61 < #include <mpi.h>
62 < #endif
64 > #include "flucq/FluctuatingChargeForces.hpp"
65  
66   namespace OpenMD {
67    
68    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
67 –                                  forceField_(NULL), info_(NULL),
69                                    haveCutoffRadius_(false),
70                                    haveDampingAlpha_(false),
71                                    haveDielectric_(false),
72 <                                  haveElectroSplines_(false)
73 <  {}
72 >                                  haveElectroSplines_(false),
73 >                                  info_(NULL), forceField_(NULL)
74 >                                  
75 >  {
76 >    flucQ_ = new FluctuatingChargeForces(info_);
77 >  }
78    
79 +  void Electrostatic::setForceField(ForceField *ff) {
80 +    forceField_ = ff;
81 +    flucQ_->setForceField(forceField_);
82 +  }
83 +
84 +  void Electrostatic::setSimulatedAtomTypes(set<AtomType*> &simtypes) {
85 +    simTypes_ = simtypes;
86 +    flucQ_->setSimulatedAtomTypes(simTypes_);
87 +  }
88 +
89    void Electrostatic::initialize() {
90      
91      Globals* simParams_ = info_->getSimParams();
# Line 249 | Line 264 | namespace OpenMD {
264      
265      RealType b0c, b1c, b2c, b3c, b4c, b5c;
266      RealType db0c_1, db0c_2, db0c_3, db0c_4, db0c_5;
267 <    RealType a2, expTerm, invArootPi;
267 >    RealType a2, expTerm, invArootPi(0.0);
268      
269      RealType r = cutoffRadius_;
270      RealType r2 = r * r;
# Line 754 | Line 769 | namespace OpenMD {
769      Tb.zero(); // Torque on site b
770      Ea.zero(); // Electric field at site a
771      Eb.zero(); // Electric field at site b
772 +    Pa = 0.0;  // Site potential at site a
773 +    Pb = 0.0;  // Site potential at site b
774      dUdCa = 0.0; // fluctuating charge force at site a
775      dUdCb = 0.0; // fluctuating charge force at site a
776      
# Line 766 | Line 783 | namespace OpenMD {
783      // Excluded potential that is still computed for fluctuating charges
784      excluded_Pot= 0.0;
785  
769 –
786      // some variables we'll need independent of electrostatic type:
787  
788      ri = 1.0 /  *(idat.rij);
# Line 829 | Line 845 | namespace OpenMD {
845        if (idat.excluded) {
846          *(idat.skippedCharge2) += C_a;
847        } else {
848 <        // only do the field if we're not excluded:
848 >        // only do the field and site potentials if we're not excluded:
849          Eb -= C_a *  pre11_ * dv01 * rhat;
850 +        Pb += C_a *  pre11_ * v01;
851        }
852      }
853      
# Line 838 | Line 855 | namespace OpenMD {
855        D_a = *(idat.dipole1);
856        rdDa = dot(rhat, D_a);
857        rxDa = cross(rhat, D_a);
858 <      if (!idat.excluded)
858 >      if (!idat.excluded) {
859          Eb -=  pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
860 +        Pb +=  pre12_ * v11 * rdDa;
861 +      }
862 +
863      }
864      
865      if (a_is_Quadrupole) {
# Line 849 | Line 869 | namespace OpenMD {
869        rQa = rhat * Q_a;
870        rdQar = dot(rhat, Qar);
871        rxQar = cross(rhat, Qar);
872 <      if (!idat.excluded)
872 >      if (!idat.excluded) {
873          Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or
874                          + rdQar * rhat * (dv22 - 2.0*v22or));
875 +        Pb += pre14_ * (v21 * trQa + v22 * rdQar);
876 +      }
877      }
878      
879      if (b_is_Charge) {
# Line 865 | Line 887 | namespace OpenMD {
887        } else {
888          // only do the field if we're not excluded:
889          Ea += C_b *  pre11_ * dv01 * rhat;
890 +        Pa += C_b *  pre11_ * v01;
891 +
892        }
893      }
894      
# Line 872 | Line 896 | namespace OpenMD {
896        D_b = *(idat.dipole2);
897        rdDb = dot(rhat, D_b);
898        rxDb = cross(rhat, D_b);
899 <      if (!idat.excluded)
899 >      if (!idat.excluded) {
900          Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b);
901 +        Pa += pre12_ * v11 * rdDb;
902 +      }
903      }
904      
905      if (b_is_Quadrupole) {
# Line 883 | Line 909 | namespace OpenMD {
909        rQb = rhat * Q_b;
910        rdQbr = dot(rhat, Qbr);
911        rxQbr = cross(rhat, Qbr);
912 <      if (!idat.excluded)
912 >      if (!idat.excluded) {
913          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
914                          + rdQbr * rhat * (dv22 - 2.0*v22or));
915 +        Pa += pre14_ * (v21 * trQb + v22 * rdQbr);
916 +      }
917      }
918 <    
918 >        
919 >
920      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
921        J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
922      }    
923 <    
923 >
924      if (a_is_Charge) {    
925        
926        if (b_is_Charge) {
927          pref =  pre11_ * *(idat.electroMult);      
928          U  += C_a * C_b * pref * v01;
929          F  += C_a * C_b * pref * dv01 * rhat;
930 <        
930 >
931          // If this is an excluded pair, there are still indirect
932          // interactions via the reaction field we must worry about:
933  
# Line 907 | Line 936 | namespace OpenMD {
936            indirect_Pot += rfContrib;
937            indirect_F   += rfContrib * 2.0 * ri * rhat;
938          }
939 <        
939 >
940          // Fluctuating charge forces are handled via Coulomb integrals
941          // for excluded pairs (i.e. those connected via bonds) and
942          // with the standard charge-charge interaction otherwise.
943  
944 <        if (idat.excluded) {          
944 >        if (idat.excluded) {
945            if (a_is_Fluctuating || b_is_Fluctuating) {
946              coulInt = J->getValueAt( *(idat.rij) );
947 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
948 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
949 <            excluded_Pot += C_a * C_b * coulInt;
921 <          }          
947 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
948 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
949 >          }
950          } else {
951            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
952 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
953 <        }
952 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
953 >        }              
954        }
955  
956        if (b_is_Dipole) {
# Line 988 | Line 1016 | namespace OpenMD {
1016          F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
1017          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
1018          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
991 –
1019          // Even if we excluded this pair from direct interactions, we
1020          // still have the reaction-field-mediated dipole-dipole
1021          // interaction:
# Line 1048 | Line 1075 | namespace OpenMD {
1075          trQaQb = QaQb.trace();
1076          rQaQb = rhat * QaQb;
1077          QaQbr = QaQb * rhat;
1078 <        QaxQb = cross(Q_a, Q_b);
1078 >        QaxQb = mCross(Q_a, Q_b);
1079          rQaQbr = dot(rQa, Qbr);
1080          rQaxQbr = cross(rQa, Qbr);
1081          
# Line 1079 | Line 1106 | namespace OpenMD {
1106          //             + 4.0 * cross(rhat, QbQar)
1107  
1108          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1082 –
1109        }
1110      }
1111  
1112      if (idat.doElectricField) {
1113        *(idat.eField1) += Ea * *(idat.electroMult);
1114        *(idat.eField2) += Eb * *(idat.electroMult);
1115 +    }
1116 +
1117 +    if (idat.doSitePotential) {
1118 +      *(idat.sPot1) += Pa * *(idat.electroMult);
1119 +      *(idat.sPot2) += Pb * *(idat.electroMult);
1120      }
1121  
1122      if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw);
# Line 1134 | Line 1165 | namespace OpenMD {
1165      bool i_is_Quadrupole = data.is_Quadrupole;
1166      bool i_is_Fluctuating = data.is_Fluctuating;
1167      RealType C_a = data.fixedCharge;  
1168 <    RealType self(0.0), preVal, DdD, trQ, trQQ;
1168 >    RealType self(0.0), preVal, DdD(0.0), trQ, trQQ;
1169  
1170      if (i_is_Dipole) {
1171        DdD = data.dipole.lengthSquare();
# Line 1142 | Line 1173 | namespace OpenMD {
1173          
1174      if (i_is_Fluctuating) {
1175        C_a += *(sdat.flucQ);
1176 <      // dVdFQ is really a force, so this is negative the derivative
1177 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1178 <      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1179 <        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1176 >
1177 >      flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ),  
1178 >                                 (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY],
1179 >                                 *(sdat.flucQfrc) );
1180 >
1181      }
1182  
1183      switch (summationMethod_) {
# Line 1248 | Line 1280 | namespace OpenMD {
1280      vector<vector<RealType> > els;
1281      vector<vector<RealType> > ems;
1282      vector<vector<RealType> > ens;
1251 –
1283      
1284      int nMax = info_->getNAtoms();
1285      
# Line 1271 | Line 1302 | namespace OpenMD {
1302      Vector3d t( 2.0 * M_PI );
1303      t.Vdiv(t, box);
1304  
1274 –    
1305      SimInfo::MoleculeIterator mi;
1306      Molecule::AtomIterator ai;
1307      int i;
# Line 1441 | Line 1471 | namespace OpenMD {
1471                    dks[i] = dk * skr[i];
1472                  }
1473                  if (data.is_Quadrupole) {
1474 <                  Q = atom->getQuadrupole();
1475 <                  Q *= mPoleConverter;
1446 <                  Qk = Q * kVec;
1474 >                  Q = atom->getQuadrupole() * mPoleConverter;
1475 >                  Qk = Q * kVec;                  
1476                    qk = dot(kVec, Qk);
1477 <                  qxk[i] = cross(kVec, Qk);
1477 >                  qxk[i] = -cross(kVec, Qk);
1478                    qkc[i] = qk * ckr[i];
1479                    qks[i] = qk * skr[i];
1480                  }              
# Line 1462 | Line 1491 | namespace OpenMD {
1491              qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1492              
1493   #ifdef IS_MPI
1494 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
1495 <                                      MPI::SUM);
1496 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE,
1497 <                                      MPI::SUM);
1498 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE,
1499 <                                      MPI::SUM);
1500 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE,
1501 <                                      MPI::SUM);
1502 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE,
1503 <                                      MPI::SUM);
1504 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE,
1505 <                                      MPI::SUM);
1494 >            MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE,
1495 >                          MPI_SUM, MPI_COMM_WORLD);
1496 >            MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE,
1497 >                          MPI_SUM, MPI_COMM_WORLD);
1498 >            MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE,
1499 >                          MPI_SUM, MPI_COMM_WORLD);
1500 >            MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE,
1501 >                          MPI_SUM, MPI_COMM_WORLD);
1502 >            MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE,
1503 >                          MPI_SUM, MPI_COMM_WORLD);
1504 >            MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE,
1505 >                          MPI_SUM, MPI_COMM_WORLD);
1506   #endif        
1507              
1508              // Accumulate potential energy and virial contribution:
# Line 1504 | Line 1533 | namespace OpenMD {
1533                  RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1534                                           -ckr[i]*(ckss+dkcs-qkss));
1535                  RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1536 <                                             +skr[i]*(ckss+dkcs-qkss));
1536 >                                            +skr[i]*(ckss+dkcs-qkss));
1537                
1538                  atom->addFrc( 4.0 * rvol * qfrc * kVec );
1539 <                
1539 >
1540 >                if (atom->isFluctuatingCharge()) {
1541 >                  atom->addFlucQFrc( - 2.0 * rvol * qtrq2 );
1542 >                }
1543 >                  
1544                  if (data.is_Dipole) {
1545                    atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1546                  }
# Line 1524 | Line 1557 | namespace OpenMD {
1557      }
1558      pot += kPot;  
1559    }
1560 +
1561 +  void Electrostatic::getSitePotentials(Atom* a1, Atom* a2, bool excluded,
1562 +                                        RealType &spot1, RealType &spot2) {
1563 +
1564 +    if (!initialized_) {
1565 +      cerr << "initializing\n";
1566 +      initialize();
1567 +      cerr << "done\n";
1568 +    }
1569 +
1570 +    const RealType mPoleConverter = 0.20819434;
1571 +
1572 +    AtomType* atype1 = a1->getAtomType();
1573 +    AtomType* atype2 = a2->getAtomType();
1574 +    int atid1 = atype1->getIdent();
1575 +    int atid2 = atype2->getIdent();
1576 +    data1 = ElectrostaticMap[Etids[atid1]];
1577 +    data2 = ElectrostaticMap[Etids[atid2]];
1578 +
1579 +    Pa = 0.0;  // Site potential at site a
1580 +    Pb = 0.0;  // Site potential at site b
1581 +
1582 +    Vector3d d = a2->getPos() - a1->getPos();
1583 +    info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(d);
1584 +    RealType rij = d.length();
1585 +    // some variables we'll need independent of electrostatic type:
1586 +
1587 +    RealType ri = 1.0 /  rij;
1588 +    rhat = d  * ri;
1589 +      
1590 +
1591 +    if ((rij >= cutoffRadius_) || excluded) {
1592 +      spot1 = 0.0;
1593 +      spot2 = 0.0;
1594 +      return;
1595 +    }
1596 +
1597 +    // logicals
1598 +
1599 +    a_is_Charge = data1.is_Charge;
1600 +    a_is_Dipole = data1.is_Dipole;
1601 +    a_is_Quadrupole = data1.is_Quadrupole;
1602 +    a_is_Fluctuating = data1.is_Fluctuating;
1603 +
1604 +    b_is_Charge = data2.is_Charge;
1605 +    b_is_Dipole = data2.is_Dipole;
1606 +    b_is_Quadrupole = data2.is_Quadrupole;
1607 +    b_is_Fluctuating = data2.is_Fluctuating;
1608 +
1609 +    // Obtain all of the required radial function values from the
1610 +    // spline structures:
1611 +    
1612 +
1613 +    if (a_is_Charge || b_is_Charge) {
1614 +      v01 = v01s->getValueAt(rij);
1615 +    }
1616 +    if (a_is_Dipole || b_is_Dipole) {
1617 +      v11 = v11s->getValueAt(rij);
1618 +      v11or = ri * v11;
1619 +    }
1620 +    if (a_is_Quadrupole || b_is_Quadrupole) {
1621 +      v21 = v21s->getValueAt(rij);
1622 +      v22 = v22s->getValueAt(rij);
1623 +      v22or = ri * v22;
1624 +    }      
1625 +
1626 +    if (a_is_Charge) {
1627 +      C_a = data1.fixedCharge;
1628 +      
1629 +      if (a_is_Fluctuating) {
1630 +        C_a += a1->getFlucQPos();
1631 +      }
1632 +      
1633 +      Pb += C_a *  pre11_ * v01;      
1634 +    }
1635 +    
1636 +    if (a_is_Dipole) {
1637 +      D_a = a1->getDipole() * mPoleConverter;
1638 +      rdDa = dot(rhat, D_a);
1639 +      Pb +=  pre12_ * v11 * rdDa;      
1640 +    }
1641 +    
1642 +    if (a_is_Quadrupole) {
1643 +      Q_a = a1->getQuadrupole() * mPoleConverter;
1644 +      trQa =  Q_a.trace();
1645 +      Qar =   Q_a * rhat;
1646 +      rdQar = dot(rhat, Qar);
1647 +      Pb += pre14_ * (v21 * trQa + v22 * rdQar);      
1648 +    }
1649 +    
1650 +    if (b_is_Charge) {
1651 +      C_b = data2.fixedCharge;
1652 +      
1653 +      if (b_is_Fluctuating)
1654 +        C_b += a2->getFlucQPos();
1655 +      
1656 +      Pa += C_b *  pre11_ * v01;
1657 +    }
1658 +    
1659 +    if (b_is_Dipole) {
1660 +      D_b = a2->getDipole() * mPoleConverter;
1661 +      rdDb = dot(rhat, D_b);
1662 +      Pa += pre12_ * v11 * rdDb;
1663 +    }
1664 +    
1665 +    if (b_is_Quadrupole) {
1666 +      Q_a = a2->getQuadrupole() * mPoleConverter;
1667 +      trQb =  Q_b.trace();
1668 +      Qbr =   Q_b * rhat;
1669 +      rdQbr = dot(rhat, Qbr);
1670 +      Pa += pre14_ * (v21 * trQb + v22 * rdQbr);      
1671 +    }
1672 +
1673 +    spot1 = Pa;
1674 +    spot2 = Pb;
1675 +  }
1676   }

Diff Legend

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