Line data    Source code

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