Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "BridgeVessel.h"
23 : #include "tools/Matrix.h"
24 : #include "core/ActionWithArguments.h"
25 : #include "StoreDataVessel.h"
26 :
27 : namespace PLMD {
28 : namespace vesselbase {
29 :
30 46 : BridgeVessel::BridgeVessel( const VesselOptions& da ):
31 : Vessel(da),
32 : inum(0),
33 : // in_normal_calculate(false)
34 : myOutputAction(NULL),
35 : myOutputValues(NULL),
36 92 : my_tmp_val(0,0)
37 : {
38 46 : }
39 :
40 513 : void BridgeVessel::resize() {
41 513 : if( myOutputAction->checkNumericalDerivatives() ) {
42 30 : mynumerical_values.resize( getAction()->getNumberOfDerivatives()*myOutputValues->getNumberOfComponents() );
43 10 : inum=0;
44 : }
45 : // This bit ensures that we can store data in a bridge function if needs be
46 : // Notice we need to resize der_list in the underlying action as this is called
47 : // from a bridge
48 513 : if( myOutputAction->mydata ) {
49 : unsigned dsize=(myOutputAction->mydata)->getSizeOfDerivativeList();
50 12 : if( getAction()->der_list.size()!=dsize ) getAction()->der_list.resize( dsize );
51 : }
52 513 : unsigned tmp=0; resizeBuffer( myOutputAction->getSizeOfBuffer( tmp ) );
53 513 : }
54 :
55 46 : void BridgeVessel::setOutputAction( ActionWithVessel* myact ) {
56 46 : ActionWithValue* checkme=dynamic_cast<ActionWithValue*>( getAction() );
57 46 : plumed_massert( checkme, "vessel in bridge must inherit from ActionWithValue");
58 :
59 46 : myOutputAction=myact;
60 46 : myOutputValues=dynamic_cast<ActionWithValue*>( myact );
61 46 : plumed_massert( myOutputValues, "bridging vessel must inherit from ActionWithValue");
62 46 : }
63 :
64 0 : std::string BridgeVessel::description() {
65 0 : plumed_merror("I shouldn't end up here");
66 : }
67 :
68 3025 : void BridgeVessel::prepare() {
69 3025 : myOutputAction->doJobsRequiredBeforeTaskList();
70 3025 : }
71 :
72 3025 : void BridgeVessel::setBufferStart( unsigned& start ) {
73 3025 : unsigned tmp=myOutputAction->getSizeOfBuffer( start );
74 3025 : }
75 :
76 102743 : MultiValue& BridgeVessel::transformDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) {
77 205022 : if( outvals.getNumberOfValues()!=myOutputAction->getNumberOfQuantities() ||
78 102279 : outvals.getNumberOfDerivatives()!=myOutputAction->getNumberOfDerivatives() ) {
79 2287 : outvals.resize( myOutputAction->getNumberOfQuantities(), myOutputAction->getNumberOfDerivatives() );
80 : }
81 102743 : myOutputAction->transformBridgedDerivatives( current, invals, outvals );
82 102743 : return outvals;
83 : }
84 :
85 102743 : void BridgeVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
86 : // in_normal_calculate=true;
87 205486 : if( myvals.get(0)<myOutputAction->getTolerance() ) return;
88 62508 : myOutputAction->calculateAllVessels( current, myvals, myvals, buffer, der_list );
89 62508 : return;
90 : }
91 :
92 3025 : void BridgeVessel::finish( const std::vector<double>& buffer ) {
93 3025 : myOutputAction->finishComputations( buffer );
94 3025 : if( myOutputAction->checkNumericalDerivatives() ) {
95 3860 : if ( inum<mynumerical_values.size() ) {
96 5400 : for(int i=0; i<myOutputValues->getNumberOfComponents(); ++i) {
97 2160 : mynumerical_values[inum]=myOutputValues->getOutputQuantity(i);
98 1080 : inum++;
99 : }
100 : plumed_dbg_assert( inum<=mynumerical_values.size() );
101 : } else {
102 850 : plumed_assert( inum==mynumerical_values.size() );
103 : }
104 : }
105 3025 : }
106 :
107 65 : void BridgeVessel::completeNumericalDerivatives() {
108 130 : unsigned nextra = myOutputAction->getNumberOfDerivatives() - getAction()->getNumberOfDerivatives();
109 130 : Matrix<double> tmpder( myOutputValues->getNumberOfComponents(), nextra );
110 65 : ActionWithVessel* vval=dynamic_cast<ActionWithVessel*>( myOutputAction );
111 1505 : for(unsigned i=0; i<nextra; ++i) {
112 1440 : vval->bridgeVariable=i; getAction()->calculate();
113 3600 : for(int j=0; j<myOutputValues->getNumberOfComponents(); ++j) tmpder(j,i) = myOutputValues->getOutputQuantity(j);
114 : }
115 130 : vval->bridgeVariable=nextra; getAction()->calculate();
116 130 : plumed_assert( inum==mynumerical_values.size() ); inum=0; // Reset inum now that we have finished calling calculate
117 130 : std::vector<double> base( myOutputValues->getNumberOfComponents() );
118 325 : for(int j=0; j<myOutputValues->getNumberOfComponents(); ++j) base[j] = myOutputValues->getOutputQuantity(j);
119 :
120 : const double delta=sqrt(epsilon);
121 65 : ActionAtomistic* aa=dynamic_cast<ActionAtomistic*>( getAction() );
122 65 : unsigned nvals=myOutputValues->getNumberOfComponents();
123 130 : for(unsigned j=0; j<nvals; ++j) ( myOutputValues->copyOutput(j) )->clearDerivatives();
124 :
125 65 : if( aa ) {
126 65 : ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( getAction() );
127 130 : plumed_assert( !aarg ); Tensor box=aa->getBox();
128 : unsigned natoms=aa->getNumberOfAtoms();
129 130 : for(unsigned j=0; j<nvals; ++j) {
130 65 : double ref=( myOutputValues->copyOutput(j) )->get();
131 65 : if( ( myOutputValues->copyOutput(j) )->getNumberOfDerivatives()>0 ) {
132 1055 : for(unsigned i=0; i<3*natoms; ++i) {
133 990 : double d=( mynumerical_values[i*nvals+j] - ref)/delta;
134 495 : ( myOutputValues->copyOutput(j) )->addDerivative(i,d);
135 : }
136 65 : Tensor virial;
137 845 : for(int i=0; i<3; i++) for(int k=0; k<3; k++) {
138 1170 : virial(i,k)=( mynumerical_values[ nvals*(3*natoms + 3*i + k) + j ]-ref)/delta;
139 : }
140 65 : virial=-matmul(box.transpose(),virial);
141 845 : for(int i=0; i<3; i++) for(int k=0; k<3; k++) ( myOutputValues->copyOutput(j) )->addDerivative(3*natoms+3*k+i,virial(k,i));
142 : }
143 : }
144 : } else {
145 0 : plumed_merror("not implemented or tested yet");
146 : // unsigned nder=myOutputAction->getNumberOfDerivatives();
147 : // for(unsigned j=0;j<nvals;++j){
148 : // double ref=( myOutputValues->copyOutput(j) )->get();
149 : // for(unsigned i=0;i<nder;++i){
150 : // double d=( mynumerical_values[i*nvals+j] - ref)/delta;
151 : // ( myOutputValues->copyOutput(j) )->addDerivative(i,d);
152 : // }
153 : // }
154 : // }
155 : }
156 : // Add the derivatives wrt to the local quantities we are working with
157 260 : for(unsigned j=0; j<nvals; ++j) {
158 : unsigned k=0;
159 785 : for(unsigned i=getAction()->getNumberOfDerivatives(); i<myOutputAction->getNumberOfDerivatives(); ++i) {
160 2880 : ( myOutputValues->copyOutput(j) )->addDerivative( i, (tmpder(j,k)-base[j])/sqrt(epsilon) ); k++;
161 : }
162 : }
163 65 : }
164 :
165 795 : bool BridgeVessel::applyForce( std::vector<double>& outforces ) {
166 1590 : bool hasforce=false; outforces.assign(outforces.size(),0.0);
167 795 : unsigned ndertot = myOutputAction->getNumberOfDerivatives();
168 795 : unsigned nextra = ndertot - getAction()->getNumberOfDerivatives();
169 795 : std::vector<double> forces( ndertot ), eforces( nextra, 0.0 );
170 4395 : for(unsigned i=0; i<myOutputAction->getNumberOfVessels(); ++i) {
171 935 : if( ( myOutputAction->getPntrToVessel(i) )->applyForce( forces ) ) {
172 : hasforce=true;
173 1043198 : for(unsigned j=0; j<outforces.size(); ++j) outforces[j]+=forces[j];
174 283 : for(unsigned j=0; j<nextra; ++j) eforces[j]+=forces[ outforces.size()+j ];
175 : }
176 : }
177 795 : if(hasforce) myOutputAction->applyBridgeForces( eforces );
178 795 : return hasforce;
179 : }
180 :
181 795 : void BridgeVessel::copyTaskFlags() {
182 795 : myOutputAction->deactivateAllTasks();
183 267018 : for(unsigned i=0; i<getAction()->nactive_tasks; ++i) myOutputAction->taskFlags[ getAction()->indexOfTaskInFullList[i] ] = 1;
184 795 : myOutputAction->lockContributors();
185 795 : }
186 :
187 13042 : MultiValue& BridgeVessel::getTemporyMultiValue() {
188 13042 : return my_tmp_val;
189 : }
190 :
191 : }
192 4839 : }
193 :
|
