Line data Source code
1 :
2 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 : Copyright (c) 2012-2023 The plumed team
4 : (see the PEOPLE file at the root of the distribution for a list of names)
5 :
6 : See http://www.plumed.org for more information.
7 :
8 : This file is part of plumed, version 2.
9 :
10 : plumed is free software: you can redistribute it and/or modify
11 : it under the terms of the GNU Lesser General Public License as published by
12 : the Free Software Foundation, either version 3 of the License, or
13 : (at your option) any later version.
14 :
15 : plumed is distributed in the hope that it will be useful,
16 : but WITHOUT ANY WARRANTY; without even the implied warranty of
17 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 : GNU Lesser General Public License for more details.
19 :
20 : You should have received a copy of the GNU Lesser General Public License
21 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
22 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
23 : #ifndef __PLUMED_colvar_Dipole_h
24 : #define __PLUMED_colvar_Dipole_h
25 : #include "Colvar.h"
26 : #include "ColvarInput.h"
27 :
28 : namespace PLMD {
29 : namespace colvar {
30 :
31 : template <typename T>
32 : class Dipole : public Colvar {
33 : std::vector<AtomNumber> ga_lista;
34 : enum /*class*/ modes:unsigned {
35 : components=1,
36 : scalar=0
37 : } mode{scalar};
38 : bool nopbc;
39 : std::vector<double> value;
40 : std::vector<double> derivs;
41 : Value* valuex=nullptr;
42 : Value* valuey=nullptr;
43 : Value* valuez=nullptr;
44 : public:
45 : using precision=T;
46 : explicit Dipole(const ActionOptions&);
47 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
48 : static unsigned getModeAndSetupValues( ActionWithValue* av );
49 : void calculate() override;
50 : static void registerKeywords(Keywords& keys);
51 : static void calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout );
52 : };
53 :
54 : template <typename T>
55 204 : void Dipole<T>::registerKeywords(Keywords& keys) {
56 204 : Colvar::registerKeywords(keys);
57 408 : keys.setDisplayName("DIPOLE");
58 204 : keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for");
59 204 : 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");
60 408 : keys.addOutputComponent("x","COMPONENTS","scalar/vector","the x-component of the dipole");
61 408 : keys.addOutputComponent("y","COMPONENTS","scalar/vector","the y-component of the dipole");
62 408 : keys.addOutputComponent("z","COMPONENTS","scalar/vector","the z-component of the dipole");
63 204 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
64 204 : keys.add("hidden","NEEDS_CHARGES","Inform the parallel task manager that we need the charges to compute this CV");
65 408 : keys.setValueDescription("scalar/vector","the DIPOLE for these atoms");
66 408 : keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
67 204 : }
68 :
69 : template <typename T>
70 60 : Dipole<T>::Dipole(const ActionOptions&ao):
71 60 : PLUMED_COLVAR_INIT(ao) {
72 60 : parseAtomList(-1,ga_lista,this);
73 60 : mode=static_cast<modes>(getModeAndSetupValues(this));
74 60 : if( mode==modes::components ) {
75 4 : value.resize(3);
76 4 : valuex=getPntrToComponent("x");
77 4 : valuey=getPntrToComponent("y");
78 4 : valuez=getPntrToComponent("z");
79 : } else {
80 56 : value.resize(1);
81 56 : derivs.resize(1,ga_lista.size());
82 : }
83 60 : parseFlag("NOPBC",nopbc);
84 60 : checkRead();
85 :
86 60 : if(nopbc) {
87 38 : log.printf(" without periodic boundary conditions\n");
88 : } else {
89 22 : log.printf(" using periodic boundary conditions\n");
90 : }
91 :
92 60 : requestAtoms(ga_lista);
93 60 : }
94 :
95 : template <typename T>
96 578 : void Dipole<T>::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
97 1156 : aa->parseAtomList("GROUP",num,t);
98 578 : if( t.size()>0 ) {
99 575 : aa->log.printf(" of %u atoms\n",static_cast<unsigned>(t.size()));
100 2489 : for(unsigned int i=0; i<t.size(); ++i) {
101 1914 : aa->log.printf(" %d", t[i].serial());
102 : }
103 575 : aa->log.printf(" \n");
104 : }
105 578 : }
106 :
107 : template <typename T>
108 63 : unsigned Dipole<T>::getModeAndSetupValues( ActionWithValue* av ) {
109 : bool c;
110 63 : av->parseFlag("COMPONENTS",c);
111 63 : if( c ) {
112 12 : av->addComponentWithDerivatives("x");
113 6 : av->componentIsNotPeriodic("x");
114 12 : av->addComponentWithDerivatives("y");
115 6 : av->componentIsNotPeriodic("y");
116 12 : av->addComponentWithDerivatives("z");
117 6 : av->componentIsNotPeriodic("z");
118 6 : return modes::components;
119 : }
120 57 : av->addValueWithDerivatives();
121 57 : av->setNotPeriodic();
122 : return modes::scalar;
123 : }
124 :
125 : // calculator
126 : template <typename T>
127 951 : void Dipole<T>::calculate() {
128 951 : if( !chargesWereSet ) {
129 0 : error("charges were not set by MD code");
130 : }
131 :
132 951 : if(!nopbc) {
133 569 : makeWhole();
134 : }
135 : const unsigned N=getNumberOfAtoms();
136 :
137 951 : auto cvout= ColvarOutput<double>::createColvarOutput(value, derivs, this);
138 951 : Dipole<double>::calculateCV( ColvarInput<double>::createColvarInput( mode, getPositions(), this ),
139 : cvout );
140 951 : if(mode==modes::scalar) {
141 6353 : for(unsigned i=0; i<N; i++) {
142 5557 : setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i));
143 : }
144 796 : setBoxDerivatives(cvout.virial[0]);
145 796 : setValue(value[0]);
146 : } else {
147 1085 : for(unsigned i=0; i<N; i++) {
148 930 : setAtomsDerivatives(valuex,i,cvout.getAtomDerivatives(0,i));
149 930 : setAtomsDerivatives(valuey,i,cvout.getAtomDerivatives(1,i));
150 930 : setAtomsDerivatives(valuez,i,cvout.getAtomDerivatives(2,i));
151 : }
152 155 : setBoxDerivatives(valuex,cvout.virial[0]);
153 155 : setBoxDerivatives(valuey,cvout.virial[1]);
154 155 : setBoxDerivatives(valuez,cvout.virial[2]);
155 155 : valuex->set(value[0]);
156 155 : valuey->set(value[1]);
157 155 : valuez->set(value[2]);
158 : }
159 951 : }
160 :
161 : template <typename T>
162 7103 : void Dipole<T>::calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout ) {
163 : using V3 = VectorTyped<T,3>;
164 7103 : const unsigned N=cvin.pos.size();
165 : T ctot=0.;
166 32106 : for(unsigned i=0; i<N; ++i) {
167 25003 : ctot += cvin.charges[i];
168 : }
169 7103 : ctot/=T(N);
170 :
171 7103 : V3 dipje;
172 32106 : for(unsigned i=0; i<N; ++i) {
173 25003 : dipje += (cvin.charges[i]-ctot)*V3(cvin.pos[i][0],cvin.pos[i][1],cvin.pos[i][2]);
174 : }
175 :
176 7103 : if( cvin.mode==modes::components ) {
177 25683 : for(unsigned i=0; i<N; i++) {
178 19386 : cvout.derivs[0][i]=(cvin.charges[i]-ctot)*V3(1.0,0.0,0.0);
179 19386 : cvout.derivs[1][i]=(cvin.charges[i]-ctot)*V3(0.0,1.0,0.0);
180 19386 : cvout.derivs[2][i]=(cvin.charges[i]-ctot)*V3(0.0,0.0,1.0);
181 : }
182 25188 : for(unsigned i=0; i<3; ++i ) {
183 18891 : cvout.values[i] = dipje[i];
184 : }
185 : } else {
186 806 : cvout.values[0] = dipje.modulo();
187 806 : T idip = 1./cvout.values[0];
188 6423 : for(unsigned i=0; i<N; i++) {
189 5617 : T dfunc=(cvin.charges[i]-ctot)*idip;
190 : cvout.derivs[0][i] = dfunc*dipje;
191 : }
192 : }
193 7103 : ColvarInput<T>::setBoxDerivativesNoPbc( cvin, cvout );
194 7103 : }
195 :
196 : }
197 : }
198 : #endif //__PLUMED_colvar_Dipole_h
|
