Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "HistogramOnGrid.h"
23 : #include "tools/KernelFunctions.h"
24 :
25 : namespace PLMD {
26 : namespace gridtools {
27 :
28 36 : void HistogramOnGrid::registerKeywords( Keywords& keys ) {
29 36 : GridVessel::registerKeywords( keys );
30 36 : keys.add("compulsory","KERNEL","the type of kernel to use");
31 36 : keys.add("compulsory","BANDWIDTH","the bandwidths");
32 36 : }
33 :
34 28 : HistogramOnGrid::HistogramOnGrid( const vesselbase::VesselOptions& da ):
35 : GridVessel(da),
36 : neigh_tot(0),
37 : addOneKernelAtATime(false),
38 : bandwidths(dimension),
39 28 : discrete(false)
40 : {
41 28 : parse("KERNEL",kerneltype);
42 28 : if( kerneltype=="discrete" || kerneltype=="DISCRETE" ) {
43 4 : discrete=true; setNoDerivatives();
44 : } else {
45 24 : parseVector("BANDWIDTH",bandwidths);
46 : }
47 28 : }
48 :
49 16 : bool HistogramOnGrid::noDiscreteKernels() const {
50 16 : return !discrete;
51 : }
52 :
53 32 : void HistogramOnGrid::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
54 : const std::vector<unsigned>& nbins, const std::vector<double>& spacing ) {
55 32 : GridVessel::setBounds( smin, smax, nbins, spacing );
56 32 : if( !discrete ) {
57 28 : std::vector<double> point(dimension,0);
58 56 : KernelFunctions kernel( point, bandwidths, kerneltype, false, 1.0, true ); neigh_tot=1;
59 56 : nneigh=kernel.getSupport( dx ); std::vector<double> support( kernel.getContinuousSupport() );
60 72 : for(unsigned i=0; i<dimension; ++i) {
61 44 : if( pbc[i] && 2*support[i]>getGridExtent(i) ) error("bandwidth is too large for periodic grid");
62 44 : neigh_tot *= (2*nneigh[i]+1);
63 28 : }
64 : }
65 32 : }
66 :
67 16002 : KernelFunctions* HistogramOnGrid::getKernelAndNeighbors( std::vector<double>& point, unsigned& num_neigh, std::vector<unsigned>& neighbors ) const {
68 16002 : if( discrete ) {
69 5 : num_neigh=1; for(unsigned i=0; i<dimension; ++i) point[i] += 0.5*dx[i];
70 5 : neighbors[0] = getIndex( point ); return NULL;
71 : } else {
72 15997 : KernelFunctions* kernel = new KernelFunctions( point, bandwidths, kerneltype, false, 1.0, true );
73 15993 : getNeighbors( kernel->getCenter(), nneigh, num_neigh, neighbors );
74 15998 : return kernel;
75 : }
76 : }
77 :
78 27019 : std::vector<Value*> HistogramOnGrid::getVectorOfValues() const {
79 27019 : std::vector<Value*> vv;
80 107201 : for(unsigned i=0; i<dimension; ++i) {
81 80188 : vv.push_back(new Value());
82 80190 : if( pbc[i] ) vv[i]->setDomain( str_min[i], str_max[i] );
83 42713 : else vv[i]->setNotPeriodic();
84 : }
85 27013 : return vv;
86 : }
87 :
88 27018 : void HistogramOnGrid::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
89 27018 : if( addOneKernelAtATime ) {
90 : plumed_dbg_assert( myvals.getNumberOfValues()==2 );
91 11046 : std::vector<double> der( dimension );
92 11061 : for(unsigned i=0; i<dimension; ++i) der[i]=myvals.getDerivative( 1, i );
93 11044 : accumulate( getAction()->getPositionInCurrentTaskList(current), myvals.get(0), myvals.get(1), der, buffer );
94 : } else {
95 : plumed_dbg_assert( myvals.getNumberOfValues()==dimension+2 );
96 15972 : std::vector<double> point( dimension ); double weight=myvals.get(0)*myvals.get( 1+dimension );
97 15971 : for(unsigned i=0; i<dimension; ++i) point[i]=myvals.get( 1+i );
98 :
99 : // Get the kernel
100 31947 : unsigned num_neigh; std::vector<unsigned> neighbors;
101 31947 : std::vector<double> der( dimension );
102 15974 : KernelFunctions* kernel=getKernelAndNeighbors( point, num_neigh, neighbors );
103 :
104 15974 : if( !kernel ) {
105 0 : plumed_dbg_assert( num_neigh==1 ); der.resize(0);
106 0 : accumulate( neighbors[0], weight, 1.0, der, buffer );
107 : } else {
108 15974 : std::vector<Value*> vv( getVectorOfValues() );
109 :
110 31947 : double newval; std::vector<double> xx( dimension );
111 15011275 : for(unsigned i=0; i<num_neigh; ++i) {
112 14995302 : unsigned ineigh=neighbors[i];
113 14995474 : if( inactive( ineigh ) ) continue ;
114 10640267 : getGridPointCoordinates( ineigh, xx );
115 10631690 : for(unsigned j=0; j<dimension; ++j) vv[j]->set(xx[j]);
116 10640227 : newval = kernel->evaluate( vv, der, true );
117 10640490 : accumulate( ineigh, weight, newval, der, buffer );
118 : }
119 15973 : delete kernel;
120 31948 : for(unsigned i=0; i<dimension; ++i) delete vv[i];
121 15974 : }
122 : }
123 27022 : }
124 :
125 11268 : void HistogramOnGrid::accumulate( const unsigned& ipoint, const double& weight, const double& dens, const std::vector<double>& der, std::vector<double>& buffer ) const {
126 11268 : buffer[bufstart+nper*ipoint] += weight*dens;
127 11272 : if( der.size()>0 ) for(unsigned j=0; j<dimension; ++j) buffer[bufstart+nper*ipoint + 1 + j] += weight*der[j];
128 11261 : }
129 :
130 51 : void HistogramOnGrid::finish( const std::vector<double>& buffer ) {
131 51 : if( addOneKernelAtATime ) {
132 11072 : for(unsigned i=0; i<getAction()->getCurrentNumberOfActiveTasks(); ++i) {
133 11048 : for(unsigned j=0; j<nper; ++j) addDataElement( nper*getAction()->getActiveTask(i)+j, buffer[bufstart+i*nper+j] );
134 : }
135 : } else {
136 27 : GridVessel::finish( buffer );
137 : }
138 51 : }
139 :
140 : }
141 2523 : }
|
