Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "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 : \plumedfile
47 : d: DIPOLE GROUP=1-10
48 : PRINT FILE=output STRIDE=5 ARG=5
49 : \endplumedfile
50 :
51 : \attention
52 : If the total charge Q of the group in non zero, then a charge Q/N will be subtracted to every atom,
53 : where N is the number of atoms. This implies that the dipole (which for a charged system depends
54 : on the position) is computed on the geometric center of the group.
55 :
56 :
57 : */
58 : //+ENDPLUMEDOC
59 :
60 141 : class Dipole : public Colvar {
61 : vector<AtomNumber> ga_lista;
62 : bool components;
63 : public:
64 : explicit Dipole(const ActionOptions&);
65 : virtual void calculate();
66 : static void registerKeywords(Keywords& keys);
67 : };
68 :
69 6499 : PLUMED_REGISTER_ACTION(Dipole,"DIPOLE")
70 :
71 48 : void Dipole::registerKeywords(Keywords& keys) {
72 48 : Colvar::registerKeywords(keys);
73 192 : keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for");
74 144 : 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");
75 192 : keys.addOutputComponent("x","COMPONENTS","the x-component of the dipole");
76 192 : keys.addOutputComponent("y","COMPONENTS","the y-component of the dipole");
77 192 : keys.addOutputComponent("z","COMPONENTS","the z-component of the dipole");
78 96 : keys.remove("NOPBC");
79 48 : }
80 :
81 47 : Dipole::Dipole(const ActionOptions&ao):
82 : PLUMED_COLVAR_INIT(ao),
83 94 : components(false)
84 : {
85 94 : parseAtomList("GROUP",ga_lista);
86 94 : parseFlag("COMPONENTS",components);
87 47 : checkRead();
88 47 : if(components) {
89 6 : addComponentWithDerivatives("x"); componentIsNotPeriodic("x");
90 6 : addComponentWithDerivatives("y"); componentIsNotPeriodic("y");
91 6 : addComponentWithDerivatives("z"); componentIsNotPeriodic("z");
92 : } else {
93 45 : addValueWithDerivatives(); setNotPeriodic();
94 : }
95 :
96 94 : log.printf(" of %u atoms\n",static_cast<unsigned>(ga_lista.size()));
97 955 : for(unsigned int i=0; i<ga_lista.size(); ++i) {
98 574 : log.printf(" %d", ga_lista[i].serial());
99 : }
100 47 : log.printf(" \n");
101 47 : requestAtoms(ga_lista);
102 47 : }
103 :
104 : // calculator
105 828 : void Dipole::calculate()
106 : {
107 : double ctot=0.;
108 : unsigned N=getNumberOfAtoms();
109 828 : vector<double> charges(N);
110 828 : Vector dipje;
111 :
112 12398 : for(unsigned i=0; i<N; ++i) {
113 11570 : charges[i]=getCharge(i);
114 5785 : ctot+=charges[i];
115 : }
116 828 : ctot/=(double)N;
117 :
118 12398 : for(unsigned i=0; i<N; ++i) {
119 11570 : charges[i]-=ctot;
120 11570 : dipje += charges[i]*getPosition(i);
121 : }
122 :
123 828 : if(!components) {
124 683 : double dipole = dipje.modulo();
125 683 : double idip = 1./dipole;
126 :
127 10513 : for(unsigned i=0; i<N; i++) {
128 9830 : double dfunc=charges[i]*idip;
129 9830 : setAtomsDerivatives(i,dfunc*dipje);
130 : }
131 683 : setBoxDerivativesNoPbc();
132 683 : setValue(dipole);
133 : } else {
134 290 : Value* valuex=getPntrToComponent("x");
135 290 : Value* valuey=getPntrToComponent("y");
136 290 : Value* valuez=getPntrToComponent("z");
137 1885 : for(unsigned i=0; i<N; i++) {
138 2610 : setAtomsDerivatives(valuex,i,charges[i]*Vector(1.0,0.0,0.0));
139 870 : setAtomsDerivatives(valuey,i,charges[i]*Vector(0.0,1.0,0.0));
140 870 : setAtomsDerivatives(valuez,i,charges[i]*Vector(0.0,0.0,1.0));
141 : }
142 145 : setBoxDerivativesNoPbc(valuex);
143 145 : setBoxDerivativesNoPbc(valuey);
144 145 : setBoxDerivativesNoPbc(valuez);
145 145 : valuex->set(dipje[0]);
146 145 : valuey->set(dipje[1]);
147 145 : valuez->set(dipje[2]);
148 : }
149 828 : }
150 :
151 : }
152 4839 : }
|
