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 :
28 : namespace PLMD {
29 : namespace crystallization {
30 :
31 0 : class VectorSum : public vesselbase::FunctionVessel {
32 : private:
33 : unsigned nder;
34 : public:
35 : static void registerKeywords( Keywords& keys );
36 : static void reserveKeyword( Keywords& keys );
37 : explicit VectorSum( const vesselbase::VesselOptions& da );
38 : std::string value_descriptor();
39 : void resize();
40 : void calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
41 : void finish( const std::vector<double>& buffer );
42 : };
43 :
44 2523 : PLUMED_REGISTER_VESSEL(VectorSum,"VSUM")
45 :
46 0 : void VectorSum::registerKeywords( Keywords& keys ) {
47 0 : vesselbase::FunctionVessel::registerKeywords(keys);
48 0 : }
49 :
50 841 : void VectorSum::reserveKeyword( Keywords& keys ) {
51 841 : keys.reserve("vessel","VSUM","calculate the norm of the sum of vectors.");
52 : keys.addOutputComponent("vsum","VSUM","the norm of sum of vectors. The output component can be refererred to elsewhere in the input "
53 841 : "file by using the label.vsum");
54 841 : }
55 :
56 0 : VectorSum::VectorSum( const vesselbase::VesselOptions& da ) :
57 : FunctionVessel(da),
58 0 : nder(0)
59 : {
60 0 : }
61 :
62 0 : std::string VectorSum::value_descriptor() {
63 0 : return "the norm of the mean vector";
64 : }
65 :
66 0 : void VectorSum::resize() {
67 0 : unsigned ncomp=getAction()->getNumberOfQuantities() - 2;
68 :
69 0 : if( getAction()->derivativesAreRequired() ) {
70 0 : nder=getAction()->getNumberOfDerivatives();
71 0 : resizeBuffer( (1+nder)*ncomp ); getFinalValue()->resizeDerivatives( nder );
72 : } else {
73 0 : nder=0; resizeBuffer(ncomp);
74 : }
75 0 : }
76 :
77 0 : void VectorSum::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
78 0 : unsigned ncomp=getAction()->getNumberOfQuantities()-2;
79 :
80 0 : double weight=myvals.get(0);
81 : plumed_dbg_assert( weight>=getTolerance() );
82 0 : buffer[bufstart] += weight;
83 0 : for(unsigned i=0; i<ncomp; ++i) buffer[bufstart + i*(1+nder)] += weight*myvals.get(2+i);
84 0 : if( !getAction()->derivativesAreRequired() ) return;
85 :
86 0 : for(unsigned i=0; i<ncomp; ++i) {
87 0 : double colvar=myvals.get(2+i);
88 0 : myvals.chainRule( 2+i, i, 1, 0, weight, bufstart, buffer );
89 0 : if( diffweight ) myvals.chainRule( 0, i, 1, 0, colvar, bufstart, buffer );
90 : }
91 0 : return;
92 : }
93 :
94 0 : void VectorSum::finish( const std::vector<double>& buffer ) {
95 0 : unsigned ncomp=getAction()->getNumberOfQuantities()-2;
96 :
97 0 : double sum=0;
98 0 : for(unsigned i=0; i<ncomp; ++i) {
99 0 : double tmp = buffer[bufstart+(nder+1)*i];
100 0 : sum+=tmp*tmp;
101 : }
102 0 : double tw = 1.0 / sqrt(sum);
103 0 : setOutputValue( sqrt(sum) );
104 0 : if( !getAction()->derivativesAreRequired() ) return;
105 :
106 0 : Value* fval=getFinalValue();
107 0 : for(unsigned icomp=0; icomp<ncomp; ++icomp) {
108 0 : double tmp = buffer[icomp*(1+nder)];
109 0 : unsigned bstart = bufstart + icomp*(nder+1) + 1;
110 0 : for(unsigned jder=0; jder<nder; ++jder) fval->addDerivative( jder, tw*tmp*buffer[bstart + jder] );
111 : }
112 : }
113 :
114 : }
115 2523 : }
|