Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2018 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 : \verbatim
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 : \endverbatim
65 : (see also \ref PRINT)
66 :
67 : */
68 : //+ENDPLUMEDOC
69 :
70 16 : class DHEnergy : public CoordinationBase {
71 : double k; // Inverse Debye screening length
72 : double constant;
73 : double epsilon;
74 :
75 : public:
76 : explicit DHEnergy(const ActionOptions&);
77 : // active methods:
78 : static void registerKeywords( Keywords& keys );
79 : virtual double pairing(double distance,double&dfunc,unsigned i,unsigned j)const;
80 : };
81 :
82 2531 : PLUMED_REGISTER_ACTION(DHEnergy,"DHENERGY")
83 :
84 9 : void DHEnergy::registerKeywords( Keywords& keys ) {
85 9 : CoordinationBase::registerKeywords(keys);
86 9 : keys.add("compulsory","I","1.0","Ionic strength (M)");
87 9 : keys.add("compulsory","TEMP","300.0","Simulation temperature (K)");
88 9 : keys.add("compulsory","EPSILON","80.0","Dielectric constant of solvent");
89 9 : }
90 :
91 : /*
92 : Global constants in SI unit used in this calculation:
93 : N_A = 6.0221412927 * 10^(23) mol^(-1) : Avogadro number
94 : q = 1.60217656535 * 10^(-19) C : proton charge
95 : e_0 = 8.854187817620 * 10^(-12) C^2/(N*m^2) : vacuum's dielectric constant
96 : k_B = 1.380648813 * 10^(-23) N*m/K : Boltzmann constant
97 : In SI unit, Debye Huckel CV is defined as:
98 : DHen = \sum_i\sum_j (q_i*q_j*q^2*N_A)/(4*pi*eps*e_0) * exp(-k*|f_ij|)/(|f_ij|)
99 : + \sum_i\sum_j (q_i*q_j*q^2*N_A)/(4*pi*epp*e_0) * (1/|r_ij| - 1/|f_ij|)
100 : = (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|))
101 : (in which |f_ij| = \sqrt(|r_ij|^2+\sigma_i*\sigma_j*exp(-|r_ij|^2/4*\sigma_i*\sigma_j)),
102 : \sigma_i and \sigma_j are the effective Born radius.)
103 : For an efficient calculation, we group constants and variables into groups:
104 : constant = (q^2*N_A)/(4*pi*e_0)
105 : tmp = 1/eps*exp(-k*|f_ij|)/(|f_ij|) + 1/epp*(1/|r_ij| - 1/|f_ij|)
106 :
107 : To speed up the loop calculation, constant can be modified as followed:
108 : constant= (q^2*N_A)/(4*pi*e_0*10^(-9))*10^(-3) (kJ/mol)
109 : = ((1.60217656535*10^(-19))^2*6.0221412927*10^(23)*10^(-3))/(4*3.14159265*8.854187817620*10^(-12)*10^(-9))
110 : = 138.935458111 (kJ/mol)
111 :
112 : */
113 :
114 8 : DHEnergy::DHEnergy(const ActionOptions&ao):
115 : Action(ao),
116 : CoordinationBase(ao),
117 : k(0.0),
118 8 : constant(0.0)
119 : {
120 : double I,T;
121 8 : parse("I",I);
122 8 : parse("TEMP",T);
123 8 : parse("EPSILON",epsilon);
124 8 : checkRead();
125 8 : if( plumed.getAtoms().usingNaturalUnits() ) error("DHENERGY cannot be used for calculations performed with natural units");
126 8 : constant=138.935458111/atoms.getUnits().getEnergy()/atoms.getUnits().getLength()*atoms.getUnits().getCharge()*atoms.getUnits().getCharge();
127 8 : k=sqrt(I/(epsilon*T))*502.903741125*atoms.getUnits().getLength();
128 8 : checkRead();
129 8 : log<<" with solvent dielectric constant "<<epsilon<<"\n";
130 8 : log<<" at temperature "<<T<<" K\n";
131 8 : log<<" at ionic strength "<<I<< "M\n";
132 8 : log<<" these parameters correspond to a screening length of "<<(1.0/k)<<"\n";
133 8 : log<<" Bibliography "<<plumed.cite("Do, Carloni, Varani and Bussi, J. Chem. Theory Comput. 9, 1720 (2013)")<<" \n";
134 8 : }
135 :
136 840 : double DHEnergy::pairing(double distance2,double&dfunc,unsigned i,unsigned j)const {
137 840 : double distance=std::sqrt(distance2);
138 840 : if(getAbsoluteIndex(i)==getAbsoluteIndex(j)) {
139 0 : dfunc=0.0;
140 0 : return 0.0;
141 : }
142 840 : double invdistance=1.0/distance;
143 840 : double tmp=exp(-k*distance)*invdistance*constant*getCharge(i)*getCharge(j)/epsilon;
144 839 : double dtmp=-(k+invdistance)*tmp;
145 839 : dfunc=dtmp*invdistance;
146 839 : return tmp;
147 : }
148 :
149 : }
150 :
151 2523 : }
|
