Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "Colvar.h"
23 : #include "ActionRegister.h"
24 :
25 : #include <string>
26 : #include <cmath>
27 :
28 : using namespace std;
29 :
30 : namespace PLMD {
31 : namespace colvar {
32 :
33 : //+PLUMEDOC COLVAR DIPOLE
34 : /*
35 : Calculate the dipole moment for a group of atoms.
36 :
37 : \warning
38 : The atoms used for \ref DIPOLE calculation should be from a whole molecule.
39 : In case the molecule is broken by the host MD code, please use \ref WHOLEMOLECULES to reconstruct it before \ref DIPOLE calculation.
40 :
41 :
42 : \par Examples
43 :
44 : The following tells plumed to calculate the dipole of the group of atoms containing
45 : the atoms from 1-10 and print it every 5 steps
46 : \verbatim
47 : d: DIPOLE GROUP=1-10
48 : PRINT FILE=output STRIDE=5 ARG=5
49 : \endverbatim
50 : (see also \ref PRINT)
51 :
52 : \attention
53 : If the total charge Q of the group in non zero, then a charge Q/N will be subtracted to every atom,
54 : where N is the number of atoms. This implies that the dipole (which for a charged system depends
55 : on the position) is computed on the geometric center of the group.
56 :
57 :
58 : */
59 : //+ENDPLUMEDOC
60 :
61 70 : class Dipole : public Colvar {
62 : vector<AtomNumber> ga_lista;
63 : bool components;
64 : public:
65 : explicit Dipole(const ActionOptions&);
66 : virtual void calculate();
67 : static void registerKeywords(Keywords& keys);
68 : };
69 :
70 2558 : PLUMED_REGISTER_ACTION(Dipole,"DIPOLE")
71 :
72 36 : void Dipole::registerKeywords(Keywords& keys) {
73 36 : Colvar::registerKeywords(keys);
74 36 : keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for");
75 36 : keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the dipole separately and store them as label.x, label.y and label.z");
76 36 : keys.addOutputComponent("x","COMPONENTS","the x-component of the dipole");
77 36 : keys.addOutputComponent("y","COMPONENTS","the y-component of the dipole");
78 36 : keys.addOutputComponent("z","COMPONENTS","the z-component of the dipole");
79 36 : keys.remove("NOPBC");
80 36 : }
81 :
82 35 : Dipole::Dipole(const ActionOptions&ao):
83 : PLUMED_COLVAR_INIT(ao),
84 35 : components(false)
85 : {
86 35 : parseAtomList("GROUP",ga_lista);
87 35 : parseFlag("COMPONENTS",components);
88 35 : checkRead();
89 35 : if(components) {
90 2 : addComponentWithDerivatives("x"); componentIsNotPeriodic("x");
91 2 : addComponentWithDerivatives("y"); componentIsNotPeriodic("y");
92 2 : addComponentWithDerivatives("z"); componentIsNotPeriodic("z");
93 : } else {
94 33 : addValueWithDerivatives(); setNotPeriodic();
95 : }
96 :
97 35 : log.printf(" of %u atoms\n",static_cast<unsigned>(ga_lista.size()));
98 250 : for(unsigned int i=0; i<ga_lista.size(); ++i) {
99 215 : log.printf(" %d", ga_lista[i].serial());
100 : }
101 35 : log.printf(" \n");
102 35 : requestAtoms(ga_lista);
103 35 : }
104 :
105 : // calculator
106 756 : void Dipole::calculate()
107 : {
108 756 : double ctot=0.;
109 756 : unsigned N=getNumberOfAtoms();
110 756 : vector<double> charges(N);
111 756 : Vector dipje;
112 :
113 6109 : for(unsigned i=0; i<N; ++i) {
114 5353 : charges[i]=getCharge(i);
115 5353 : ctot+=charges[i];
116 : }
117 756 : ctot/=(double)N;
118 :
119 6109 : for(unsigned i=0; i<N; ++i) {
120 5353 : charges[i]-=ctot;
121 5353 : dipje += charges[i]*getPosition(i);
122 : }
123 :
124 756 : if(!components) {
125 611 : double dipole = dipje.modulo();
126 611 : double idip = 1./dipole;
127 :
128 5094 : for(unsigned i=0; i<N; i++) {
129 4483 : double dfunc=charges[i]*idip;
130 4483 : setAtomsDerivatives(i,dfunc*dipje);
131 : }
132 611 : setBoxDerivativesNoPbc();
133 611 : setValue(dipole);
134 : } else {
135 145 : Value* valuex=getPntrToComponent("x");
136 145 : Value* valuey=getPntrToComponent("y");
137 145 : Value* valuez=getPntrToComponent("z");
138 1015 : for(unsigned i=0; i<N; i++) {
139 870 : setAtomsDerivatives(valuex,i,charges[i]*Vector(1.0,0.0,0.0));
140 870 : setAtomsDerivatives(valuey,i,charges[i]*Vector(0.0,1.0,0.0));
141 870 : setAtomsDerivatives(valuez,i,charges[i]*Vector(0.0,0.0,1.0));
142 : }
143 145 : setBoxDerivativesNoPbc(valuex);
144 145 : setBoxDerivativesNoPbc(valuey);
145 145 : setBoxDerivativesNoPbc(valuez);
146 145 : valuex->set(dipje[0]);
147 145 : valuey->set(dipje[1]);
148 145 : valuez->set(dipje[2]);
149 756 : }
150 756 : }
151 :
152 : }
153 2523 : }
|
