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 1993 by gezelter, Tue Apr 29 17:32:31 2014 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    
# Line 69 | Line 71 | namespace OpenMD {
71                                    haveDampingAlpha_(false),
72                                    haveDielectric_(false),
73                                    haveElectroSplines_(false)
74 <  {}
74 >  {
75 >    flucQ_ = new FluctuatingChargeForces(info_);
76 >  }
77    
78 +  void Electrostatic::setForceField(ForceField *ff) {
79 +    forceField_ = ff;
80 +    flucQ_->setForceField(forceField_);
81 +  }
82 +
83 +  void Electrostatic::setSimulatedAtomTypes(set<AtomType*> &simtypes) {
84 +    simTypes_ = simtypes;
85 +    flucQ_->setSimulatedAtomTypes(simTypes_);
86 +  }
87 +
88    void Electrostatic::initialize() {
89      
90      Globals* simParams_ = info_->getSimParams();
# Line 754 | Line 768 | namespace OpenMD {
768      Tb.zero(); // Torque on site b
769      Ea.zero(); // Electric field at site a
770      Eb.zero(); // Electric field at site b
771 +    Pa = 0.0;  // Site potential at site a
772 +    Pb = 0.0;  // Site potential at site b
773      dUdCa = 0.0; // fluctuating charge force at site a
774      dUdCb = 0.0; // fluctuating charge force at site a
775      
# Line 766 | Line 782 | namespace OpenMD {
782      // Excluded potential that is still computed for fluctuating charges
783      excluded_Pot= 0.0;
784  
769 –
785      // some variables we'll need independent of electrostatic type:
786  
787      ri = 1.0 /  *(idat.rij);
# Line 829 | Line 844 | namespace OpenMD {
844        if (idat.excluded) {
845          *(idat.skippedCharge2) += C_a;
846        } else {
847 <        // only do the field if we're not excluded:
847 >        // only do the field and site potentials if we're not excluded:
848          Eb -= C_a *  pre11_ * dv01 * rhat;
849 +        Pb += C_a *  pre11_ * v01;
850        }
851      }
852      
# Line 838 | Line 854 | namespace OpenMD {
854        D_a = *(idat.dipole1);
855        rdDa = dot(rhat, D_a);
856        rxDa = cross(rhat, D_a);
857 <      if (!idat.excluded)
857 >      if (!idat.excluded) {
858          Eb -=  pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
859 +        Pb +=  pre12_ * v11 * rdDa;
860 +      }
861 +
862      }
863      
864      if (a_is_Quadrupole) {
# Line 849 | Line 868 | namespace OpenMD {
868        rQa = rhat * Q_a;
869        rdQar = dot(rhat, Qar);
870        rxQar = cross(rhat, Qar);
871 <      if (!idat.excluded)
871 >      if (!idat.excluded) {
872          Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or
873                          + rdQar * rhat * (dv22 - 2.0*v22or));
874 +        Pb += pre14_ * (v21 * trQa + v22 * rdQar);
875 +      }
876      }
877      
878      if (b_is_Charge) {
# Line 865 | Line 886 | namespace OpenMD {
886        } else {
887          // only do the field if we're not excluded:
888          Ea += C_b *  pre11_ * dv01 * rhat;
889 +        Pa += C_b *  pre11_ * v01;
890 +
891        }
892      }
893      
# Line 872 | Line 895 | namespace OpenMD {
895        D_b = *(idat.dipole2);
896        rdDb = dot(rhat, D_b);
897        rxDb = cross(rhat, D_b);
898 <      if (!idat.excluded)
898 >      if (!idat.excluded) {
899          Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b);
900 +        Pa += pre12_ * v11 * rdDb;
901 +      }
902      }
903      
904      if (b_is_Quadrupole) {
# Line 883 | Line 908 | namespace OpenMD {
908        rQb = rhat * Q_b;
909        rdQbr = dot(rhat, Qbr);
910        rxQbr = cross(rhat, Qbr);
911 <      if (!idat.excluded)
911 >      if (!idat.excluded) {
912          Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
913                          + rdQbr * rhat * (dv22 - 2.0*v22or));
914 +        Pa += pre14_ * (v21 * trQb + v22 * rdQbr);
915 +      }
916      }
917 <    
917 >        
918 >
919      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
920        J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
921      }    
922 <    
922 >
923      if (a_is_Charge) {    
924        
925        if (b_is_Charge) {
926          pref =  pre11_ * *(idat.electroMult);      
927          U  += C_a * C_b * pref * v01;
928          F  += C_a * C_b * pref * dv01 * rhat;
929 <        
929 >
930          // If this is an excluded pair, there are still indirect
931          // interactions via the reaction field we must worry about:
932  
# Line 907 | Line 935 | namespace OpenMD {
935            indirect_Pot += rfContrib;
936            indirect_F   += rfContrib * 2.0 * ri * rhat;
937          }
938 <        
938 >
939          // Fluctuating charge forces are handled via Coulomb integrals
940          // for excluded pairs (i.e. those connected via bonds) and
941          // with the standard charge-charge interaction otherwise.
942  
943 <        if (idat.excluded) {          
943 >        if (idat.excluded) {
944            if (a_is_Fluctuating || b_is_Fluctuating) {
945              coulInt = J->getValueAt( *(idat.rij) );
946 <            if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
947 <            if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
948 <            excluded_Pot += C_a * C_b * coulInt;
921 <          }          
946 >            if (a_is_Fluctuating) dUdCa += C_b * coulInt;
947 >            if (b_is_Fluctuating) dUdCb += C_a * coulInt;          
948 >          }
949          } else {
950            if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
951 <          if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
952 <        }
951 >          if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
952 >        }              
953        }
954  
955        if (b_is_Dipole) {
# Line 988 | Line 1015 | namespace OpenMD {
1015          F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
1016          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
1017          Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
991 –
1018          // Even if we excluded this pair from direct interactions, we
1019          // still have the reaction-field-mediated dipole-dipole
1020          // interaction:
# Line 1048 | Line 1074 | namespace OpenMD {
1074          trQaQb = QaQb.trace();
1075          rQaQb = rhat * QaQb;
1076          QaQbr = QaQb * rhat;
1077 <        QaxQb = cross(Q_a, Q_b);
1077 >        QaxQb = mCross(Q_a, Q_b);
1078          rQaQbr = dot(rQa, Qbr);
1079          rQaxQbr = cross(rQa, Qbr);
1080          
# Line 1079 | Line 1105 | namespace OpenMD {
1105          //             + 4.0 * cross(rhat, QbQar)
1106  
1107          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1082 –
1108        }
1109      }
1110  
# Line 1088 | Line 1113 | namespace OpenMD {
1113        *(idat.eField2) += Eb * *(idat.electroMult);
1114      }
1115  
1116 +    if (idat.doSitePotential) {
1117 +      *(idat.sPot1) += Pa * *(idat.electroMult);
1118 +      *(idat.sPot2) += Pb * *(idat.electroMult);
1119 +    }
1120 +
1121      if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw);
1122      if (b_is_Fluctuating) *(idat.dVdFQ2) += dUdCb * *(idat.sw);
1123  
# Line 1142 | Line 1172 | namespace OpenMD {
1172          
1173      if (i_is_Fluctuating) {
1174        C_a += *(sdat.flucQ);
1175 <      // dVdFQ is really a force, so this is negative the derivative
1176 <      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1177 <      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1178 <        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1175 >
1176 >      flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ),  
1177 >                                 (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY],
1178 >                                 *(sdat.flucQfrc) );
1179 >
1180      }
1181  
1182      switch (summationMethod_) {
# Line 1248 | Line 1279 | namespace OpenMD {
1279      vector<vector<RealType> > els;
1280      vector<vector<RealType> > ems;
1281      vector<vector<RealType> > ens;
1251 –
1282      
1283      int nMax = info_->getNAtoms();
1284      
# Line 1271 | Line 1301 | namespace OpenMD {
1301      Vector3d t( 2.0 * M_PI );
1302      t.Vdiv(t, box);
1303  
1274 –    
1304      SimInfo::MoleculeIterator mi;
1305      Molecule::AtomIterator ai;
1306      int i;
# Line 1441 | Line 1470 | namespace OpenMD {
1470                    dks[i] = dk * skr[i];
1471                  }
1472                  if (data.is_Quadrupole) {
1473 <                  Q = atom->getQuadrupole();
1474 <                  Q *= mPoleConverter;
1446 <                  Qk = Q * kVec;
1473 >                  Q = atom->getQuadrupole() * mPoleConverter;
1474 >                  Qk = Q * kVec;                  
1475                    qk = dot(kVec, Qk);
1476 <                  qxk[i] = cross(kVec, Qk);
1476 >                  qxk[i] = -cross(kVec, Qk);
1477                    qkc[i] = qk * ckr[i];
1478                    qks[i] = qk * skr[i];
1479                  }              
# Line 1462 | Line 1490 | namespace OpenMD {
1490              qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1491              
1492   #ifdef IS_MPI
1493 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE,
1494 <                                      MPI::SUM);
1495 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE,
1496 <                                      MPI::SUM);
1497 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE,
1498 <                                      MPI::SUM);
1499 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE,
1500 <                                      MPI::SUM);
1501 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE,
1502 <                                      MPI::SUM);
1503 <            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE,
1504 <                                      MPI::SUM);
1493 >            MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE,
1494 >                          MPI_SUM, MPI_COMM_WORLD);
1495 >            MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE,
1496 >                          MPI_SUM, MPI_COMM_WORLD);
1497 >            MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE,
1498 >                          MPI_SUM, MPI_COMM_WORLD);
1499 >            MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE,
1500 >                          MPI_SUM, MPI_COMM_WORLD);
1501 >            MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE,
1502 >                          MPI_SUM, MPI_COMM_WORLD);
1503 >            MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE,
1504 >                          MPI_SUM, MPI_COMM_WORLD);
1505   #endif        
1506              
1507              // Accumulate potential energy and virial contribution:
# Line 1504 | Line 1532 | namespace OpenMD {
1532                  RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1533                                           -ckr[i]*(ckss+dkcs-qkss));
1534                  RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1535 <                                             +skr[i]*(ckss+dkcs-qkss));
1535 >                                            +skr[i]*(ckss+dkcs-qkss));
1536                
1537                  atom->addFrc( 4.0 * rvol * qfrc * kVec );
1538 <                
1538 >
1539 >                if (atom->isFluctuatingCharge()) {
1540 >                  atom->addFlucQFrc( - 2.0 * rvol * qtrq2 );
1541 >                }
1542 >                  
1543                  if (data.is_Dipole) {
1544                    atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1545                  }

Diff Legend

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