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