Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "CoordinationBase.h"
      23             : #include "tools/SwitchingFunction.h"
      24             : #include "core/ActionRegister.h"
      25             : #include "core/PlumedMain.h"
      26             : 
      27             : namespace PLMD {
      28             : namespace colvar {
      29             : 
      30             : //+PLUMEDOC COLVAR DHENERGY
      31             : /*
      32             : Calculate Debye-Huckel interaction energy among GROUPA and GROUPB.
      33             : 
      34             : This variable calculates the electrostatic interaction among GROUPA and GROUPB
      35             : using a Debye-Huckel approximation defined as
      36             : \f[
      37             : \frac{1}{4\pi\epsilon_r\epsilon_0}
      38             : \sum_{i\in A} \sum_{j \in B} q_i q_j
      39             : \frac{e^{-\kappa |{\bf r}_{ij}|}}{|{\bf r}_{ij}|}
      40             : \f]
      41             : 
      42             : This collective variable can be used to analyze or induce electrostatically driven reactions \cite do13jctc.
      43             : Notice that the value of the DHENERGY is returned in plumed units (see \ref UNITS).
      44             : 
      45             : If GROUPB is empty, it will sum the N*(N-1)/2 pairs in GROUPA. This avoids computing
      46             : twice permuted indexes (e.g. pair (i,j) and (j,i)) thus running at twice the speed.
      47             : 
      48             : Notice that if there are common atoms between GROUPA and GROUPB their interaction is discarded.
      49             : 
      50             : 
      51             : \par Examples
      52             : 
      53             : \plumedfile
      54             : # this is printing the electrostatic interaction between two groups of atoms
      55             : dh: DHENERGY GROUPA=1-10 GROUPB=11-20 EPSILON=80.0 I=0.1 TEMP=300.0
      56             : PRINT ARG=dh
      57             : \endplumedfile
      58             : 
      59             : */
      60             : //+ENDPLUMEDOC
      61             : 
      62             : class DHEnergy : public CoordinationBase {
      63             :   double k; // Inverse Debye screening length
      64             :   double constant;
      65             :   double epsilon;
      66             : 
      67             : public:
      68             :   explicit DHEnergy(const ActionOptions&);
      69             : // active methods:
      70             :   static void registerKeywords( Keywords& keys );
      71             :   double pairing(double distance,double&dfunc,unsigned i,unsigned j)const override;
      72             : };
      73             : 
      74             : PLUMED_REGISTER_ACTION(DHEnergy,"DHENERGY")
      75             : 
      76          10 : void DHEnergy::registerKeywords( Keywords& keys ) {
      77          10 :   CoordinationBase::registerKeywords(keys);
      78          20 :   keys.add("compulsory","I","1.0","Ionic strength (M)");
      79          20 :   keys.add("compulsory","TEMP","300.0","Simulation temperature (K)");
      80          20 :   keys.add("compulsory","EPSILON","80.0","Dielectric constant of solvent");
      81          10 :   keys.setValueDescription("the value of the DHENERGY");
      82          10 : }
      83             : 
      84             : /*
      85             : Global constants in SI unit used in this calculation:
      86             :       N_A = 6.0221412927 * 10^(23) mol^(-1) : Avogadro number
      87             :       q = 1.60217656535 * 10^(-19) C : proton charge
      88             :       e_0 = 8.854187817620 * 10^(-12) C^2/(N*m^2) : vacuum's dielectric constant
      89             :       k_B = 1.380648813 * 10^(-23) N*m/K : Boltzmann constant
      90             : In SI unit, Debye Huckel CV is defined as:
      91             :       DHen = \sum_i\sum_j (q_i*q_j*q^2*N_A)/(4*pi*eps*e_0) * exp(-k*|f_ij|)/(|f_ij|)
      92             :              + \sum_i\sum_j (q_i*q_j*q^2*N_A)/(4*pi*epp*e_0) * (1/|r_ij| - 1/|f_ij|)
      93             :            = (q^2*N_A)/(4*pi*e_0) * \sum_i\sum_j q_i*q_j * (exp(-k*|f_ij|)/(eps*|f_ij|) + 1/epp*(1/|r_ij| - 1/|f_ij|))
      94             : (in which |f_ij| = \sqrt(|r_ij|^2+\sigma_i*\sigma_j*exp(-|r_ij|^2/4*\sigma_i*\sigma_j)),
      95             :  \sigma_i and \sigma_j are the effective Born radius.)
      96             : For an efficient calculation, we group constants and variables into groups:
      97             :       constant = (q^2*N_A)/(4*pi*e_0)
      98             :       tmp = 1/eps*exp(-k*|f_ij|)/(|f_ij|) + 1/epp*(1/|r_ij| - 1/|f_ij|)
      99             : 
     100             : To speed up the loop calculation, constant can be modified as followed:
     101             :       constant= (q^2*N_A)/(4*pi*e_0*10^(-9))*10^(-3) (kJ/mol)
     102             :               = ((1.60217656535*10^(-19))^2*6.0221412927*10^(23)*10^(-3))/(4*3.14159265*8.854187817620*10^(-12)*10^(-9))
     103             :               = 138.935458111 (kJ/mol)
     104             : 
     105             : */
     106             : 
     107           8 : DHEnergy::DHEnergy(const ActionOptions&ao):
     108             :   Action(ao),
     109             :   CoordinationBase(ao),
     110           8 :   k(0.0),
     111           8 :   constant(0.0)
     112             : {
     113             :   double I,T;
     114           8 :   parse("I",I);
     115           8 :   parse("TEMP",T);
     116           8 :   parse("EPSILON",epsilon);
     117           8 :   checkRead();
     118           8 :   if( usingNaturalUnits() ) error("DHENERGY cannot be used for calculations performed with natural units");
     119           8 :   constant=138.935458111/getUnits().getEnergy()/getUnits().getLength()*getUnits().getCharge()*getUnits().getCharge();
     120           8 :   k=std::sqrt(I/(epsilon*T))*502.903741125*getUnits().getLength();
     121           8 :   checkRead();
     122           8 :   log<<"  with solvent dielectric constant "<<epsilon<<"\n";
     123           8 :   log<<"  at temperature "<<T<<" K\n";
     124           8 :   log<<"  at ionic strength "<<I<< "M\n";
     125           8 :   log<<"  these parameters correspond to a screening length of "<<(1.0/k)<<"\n";
     126          16 :   log<<"  Bibliography "<<plumed.cite("Do, Carloni, Varani and Bussi, J. Chem. Theory Comput. 9, 1720 (2013)")<<" \n";
     127           8 : }
     128             : 
     129         840 : double DHEnergy::pairing(double distance2,double&dfunc,unsigned i,unsigned j)const {
     130         840 :   double distance=std::sqrt(distance2);
     131         840 :   if(getAbsoluteIndex(i)==getAbsoluteIndex(j)) {
     132           0 :     dfunc=0.0;
     133           0 :     return 0.0;
     134             :   }
     135         840 :   double invdistance=1.0/distance;
     136         840 :   double tmp=std::exp(-k*distance)*invdistance*constant*getCharge(i)*getCharge(j)/epsilon;
     137         840 :   double dtmp=-(k+invdistance)*tmp;
     138         840 :   dfunc=dtmp*invdistance;
     139         840 :   return tmp;
     140             : }
     141             : 
     142             : }
     143             : 
     144             : }