Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 "vesselbase/VesselRegister.h"
23 : #include "vesselbase/FunctionVessel.h"
24 : #include "vesselbase/ActionWithVessel.h"
25 : #include "multicolvar/ActionVolume.h"
26 : #include "VectorMultiColvar.h"
27 : #include "Gradient.h"
28 :
29 : namespace PLMD {
30 : namespace crystallization {
31 :
32 4 : class GradientVessel : public vesselbase::FunctionVessel {
33 : private:
34 : bool isdens;
35 : unsigned nweights, ncomponents;
36 : std::vector<unsigned> starts;
37 : public:
38 : static void registerKeywords( Keywords& keys );
39 : static void reserveKeyword( Keywords& keys );
40 : explicit GradientVessel( const vesselbase::VesselOptions& da );
41 : std::string value_descriptor();
42 : void resize();
43 : void calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
44 : void finish( const std::vector<double>& buffer );
45 : };
46 :
47 2525 : PLUMED_REGISTER_VESSEL(GradientVessel,"GRADIENT")
48 :
49 2 : void GradientVessel::registerKeywords( Keywords& keys ) {
50 2 : vesselbase::FunctionVessel::registerKeywords(keys);
51 2 : }
52 :
53 841 : void GradientVessel::reserveKeyword( Keywords& keys ) {
54 841 : keys.reserve("vessel","GRADIENT","calculate the gradient");
55 841 : keys.addOutputComponent("gradient","GRADIENT","the gradient");
56 841 : }
57 :
58 2 : GradientVessel::GradientVessel( const vesselbase::VesselOptions& da ) :
59 2 : FunctionVessel(da)
60 : {
61 2 : Gradient* vg=dynamic_cast<Gradient*>( getAction() );
62 2 : plumed_assert( vg ); isdens=(vg->getPntrToMultiColvar())->isDensity();
63 2 : nweights = vg->nbins[0] + vg->nbins[1] + vg->nbins[2];
64 2 : if( (vg->getPntrToMultiColvar())->getNumberOfQuantities()>2 ) {
65 0 : ncomponents = (vg->getPntrToMultiColvar())->getNumberOfQuantities() - 2;
66 : } else {
67 2 : ncomponents = 1;
68 : }
69 :
70 2 : starts.push_back(0);
71 2 : if( vg->nbins[0]>0 ) {
72 2 : starts.push_back( vg->nbins[0] );
73 2 : if( vg->nbins[1]>0 ) {
74 0 : starts.push_back( vg->nbins[0] + vg->nbins[1] );
75 0 : if( vg->nbins[2]>0 ) starts.push_back( nweights );
76 2 : } else if( vg->nbins[2]>0 ) {
77 0 : starts.push_back( nweights );
78 : }
79 0 : } else if( vg->nbins[1]>0 ) {
80 0 : starts.push_back( vg->nbins[1] );
81 0 : if( vg->nbins[2]>0 ) starts.push_back( nweights );
82 0 : } else if( vg->nbins[2]>0 ) {
83 0 : starts.push_back( nweights );
84 : }
85 2 : }
86 :
87 2 : std::string GradientVessel::value_descriptor() {
88 2 : return "the gradient";
89 : }
90 :
91 4 : void GradientVessel::resize() {
92 4 : if( getAction()->derivativesAreRequired() ) {
93 2 : unsigned nder=getAction()->getNumberOfDerivatives();
94 2 : resizeBuffer( (1+nder)*(ncomponents+1)*nweights );
95 2 : getFinalValue()->resizeDerivatives( nder );
96 : } else {
97 2 : resizeBuffer( (ncomponents+1)*nweights );
98 : }
99 4 : }
100 :
101 11570 : void GradientVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
102 : unsigned nder;
103 11570 : if( getAction()->derivativesAreRequired() ) nder=getAction()->getNumberOfDerivatives();
104 0 : else nder=0;
105 11574 : unsigned wstart, cstart; if( ncomponents==1 ) { cstart=1; wstart=2; } else { cstart=2; wstart=2+ncomponents; }
106 :
107 57972 : for(unsigned iw=0; iw<nweights; ++iw) {
108 46378 : unsigned xx = (ncomponents+1)*iw;
109 46378 : double weight=myvals.get(wstart+iw);
110 46398 : buffer[bufstart+xx*(nder+1)] += weight;
111 46398 : myvals.chainRule( wstart + iw, xx, 1, 0, 1.0, bufstart, buffer );
112 92778 : for(unsigned jc=0; jc<ncomponents; ++jc) {
113 46380 : double colvar=myvals.get( cstart + jc );
114 46381 : buffer[bufstart+(xx+1+jc)*(nder+1) ] += weight*colvar;
115 46386 : myvals.chainRule( cstart + jc, xx + 1 + jc, 1, 0, weight, bufstart, buffer );
116 46394 : myvals.chainRule( wstart + iw, xx + 1 + jc, 1, 0, colvar, bufstart, buffer );
117 : }
118 : }
119 11594 : }
120 :
121 232 : void GradientVessel::finish( const std::vector<double>& buffer ) {
122 232 : std::vector<double> val_interm( ncomponents*nweights );
123 : unsigned nder;
124 232 : if( getAction()->derivativesAreRequired() ) nder=getAction()->getNumberOfDerivatives();
125 0 : else nder=0;
126 464 : Matrix<double> der_interm( ncomponents*nweights, nder ); der_interm=0;
127 :
128 232 : if( isdens ) {
129 580 : for(unsigned iw=0; iw<nweights; ++iw) {
130 464 : val_interm[iw] = buffer[bufstart + 2*iw*(1+nder)];
131 464 : if( getAction()->derivativesAreRequired() ) {
132 464 : unsigned wstart = bufstart + 2*iw*(nder+1) + 1;
133 464 : for(unsigned jder=0; jder<nder; ++jder) der_interm( iw, jder ) += buffer[ wstart + jder ];
134 : }
135 : }
136 : } else {
137 580 : for(unsigned iw=0; iw<nweights; ++iw) {
138 464 : unsigned xx = (ncomponents+1)*iw;
139 464 : double ww=buffer[bufstart + xx*(1+nder)];
140 464 : for(unsigned jc=0; jc<ncomponents; ++jc) val_interm[ iw*ncomponents + jc ] = buffer[bufstart + (xx+1+jc)*(1+nder)] / ww;
141 464 : if( getAction()->derivativesAreRequired() ) {
142 464 : unsigned wstart = bufstart + xx*(nder+1) + 1;
143 928 : for(unsigned jc=0; jc<ncomponents; ++jc) {
144 464 : unsigned bstart = bufstart + ( xx + 1 + jc )*(nder+1) + 1;
145 464 : double val = buffer[bufstart + (nder+1)*(xx+1+jc)];
146 75632 : for(unsigned jder=0; jder<nder; ++jder)
147 75168 : der_interm( iw*ncomponents + jc, jder ) = (1.0/ww)*buffer[bstart + jder] - (val/(ww*ww))*buffer[wstart + jder];
148 : }
149 : }
150 : }
151 : }
152 :
153 232 : double tmp, diff2=0.0;
154 :
155 232 : if( getAction()->derivativesAreRequired() ) {
156 232 : Value* fval=getFinalValue();
157 464 : for(unsigned j=0; j<starts.size()-1; ++j) {
158 1160 : for(unsigned bin=starts[j]; bin<starts[j+1]; ++bin) {
159 1856 : for(unsigned jc=0; jc<ncomponents; ++jc) {
160 928 : if( bin==starts[j] ) {
161 232 : tmp=val_interm[(starts[j+1]-1)*ncomponents + jc] - val_interm[bin*ncomponents + jc];
162 37816 : for(unsigned jder=0; jder<nder; ++jder) {
163 37584 : fval->addDerivative( jder, +2.0*tmp*der_interm( (starts[j+1]-1)*ncomponents + jc, jder) );
164 37584 : fval->addDerivative( jder, -2.0*tmp*der_interm( bin*ncomponents + jc, jder ) );
165 : }
166 : } else {
167 696 : tmp=val_interm[(bin-1)*ncomponents + jc] - val_interm[bin*ncomponents + jc];
168 113448 : for(unsigned jder=0; jder<nder; ++jder) {
169 112752 : fval->addDerivative( jder, +2.0*tmp*der_interm( (bin-1)*ncomponents + jc, jder) );
170 112752 : fval->addDerivative( jder, -2.0*tmp*der_interm( bin*ncomponents + jc, jder ) );
171 : }
172 : }
173 928 : diff2+=tmp*tmp;
174 : }
175 : }
176 : }
177 : } else {
178 0 : for(unsigned j=0; j<starts.size()-1; ++j) {
179 0 : for(unsigned bin=starts[j]; bin<starts[j+1]; ++bin) {
180 0 : for(unsigned jc=0; jc<ncomponents; ++jc) {
181 0 : if( bin==starts[j] ) tmp=val_interm[(starts[j+1]-1)*ncomponents + jc] - val_interm[bin*ncomponents + jc];
182 0 : else tmp=val_interm[(bin-1)*ncomponents + jc] - val_interm[bin*ncomponents + jc];
183 0 : diff2+=tmp*tmp;
184 : }
185 : }
186 : }
187 : }
188 464 : setOutputValue( diff2 );
189 232 : }
190 :
191 : }
192 2523 : }
|
