Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2019 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 6460 : 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 4839 : }
|
