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 "ReferenceArguments.h"
23 : #include "ReferenceAtoms.h"
24 : #include "tools/OFile.h"
25 : #include "core/Value.h"
26 :
27 : namespace PLMD {
28 :
29 1162 : ReferenceArguments::ReferenceArguments( const ReferenceConfigurationOptions& ro ):
30 : ReferenceConfiguration(ro),
31 : hasweights(false),
32 2324 : hasmetric(false)
33 : {
34 1162 : }
35 :
36 100 : void ReferenceArguments::readArgumentsFromPDB( const PDB& pdb ) {
37 100 : ReferenceAtoms* aref=dynamic_cast<ReferenceAtoms*>( this );
38 192 : if( !aref ) parseVector( "ARG", arg_names );
39 16 : else parseVector( "ARG", arg_names, true );
40 :
41 200 : reference_args.resize( arg_names.size() );
42 752 : for(unsigned i=0; i<arg_names.size(); ++i) parse( arg_names[i], reference_args[i] );
43 :
44 100 : if( hasweights ) {
45 0 : plumed_massert( !hasmetric, "should not have weights if we are using metric");
46 0 : weights.resize( arg_names.size() ); sqrtweight.resize( arg_names.size() );
47 0 : for(unsigned i=0; i<reference_args.size(); ++i) {
48 0 : parse( "sigma_" + arg_names[i], weights[i] );
49 0 : sqrtweight[i] = sqrt( weights[i] );
50 : }
51 100 : } else if( hasmetric ) {
52 : plumed_massert( !hasweights, "should not have weights if we are using metric");
53 0 : double thissig; metric.resize( arg_names.size(), arg_names.size() );
54 0 : for(unsigned i=0; i<reference_args.size(); ++i) {
55 0 : for(unsigned j=i; j<reference_args.size(); ++j) {
56 0 : parse( "sigma_" + arg_names[i] + "_" + arg_names[j], thissig );
57 0 : metric(i,j)=metric(j,i)=thissig;
58 : }
59 : }
60 : } else {
61 200 : weights.resize( arg_names.size() ); sqrtweight.resize( arg_names.size() );
62 936 : for(unsigned i=0; i<weights.size(); ++i) sqrtweight[i]=weights[i]=1.0;
63 : }
64 100 : }
65 :
66 1060 : void ReferenceArguments::setArgumentNames( const std::vector<std::string>& arg_vals ) {
67 2120 : reference_args.resize( arg_vals.size() );
68 2120 : arg_names.resize( arg_vals.size() );
69 2120 : arg_der_index.resize( arg_vals.size() );
70 8222 : for(unsigned i=0; i<arg_vals.size(); ++i) {
71 2034 : arg_names[i]=arg_vals[i]; arg_der_index[i]=i;
72 : }
73 1060 : if( hasmetric ) metric.resize( arg_vals.size(), arg_vals.size() );
74 1060 : else weights.resize( arg_vals.size() );
75 1060 : }
76 :
77 1922 : void ReferenceArguments::setReferenceArguments( const std::vector<double>& arg_vals, const std::vector<double>& sigma ) {
78 1922 : moveReferenceArguments( arg_vals );
79 :
80 1922 : if( hasmetric ) {
81 : unsigned k=0;
82 0 : for(unsigned i=0; i<reference_args.size(); ++i) {
83 0 : for(unsigned j=i; j<reference_args.size(); ++j) {
84 0 : metric(i,j)=metric(j,i)=sigma[k]; k++;
85 : }
86 : }
87 0 : plumed_assert( k==sigma.size() );
88 : } else {
89 1922 : plumed_assert( reference_args.size()==sigma.size() );
90 15262 : for(unsigned i=0; i<reference_args.size(); ++i) weights[i]=sigma[i];
91 : }
92 1922 : }
93 :
94 1922 : void ReferenceArguments::moveReferenceArguments( const std::vector<double>& arg_vals ) {
95 : plumed_dbg_assert( reference_args.size()==arg_vals.size() );
96 15262 : for(unsigned i=0; i<arg_vals.size(); ++i) reference_args[i]=arg_vals[i];
97 1922 : }
98 :
99 100 : void ReferenceArguments::getArgumentRequests( std::vector<std::string>& argout, bool disable_checks ) {
100 200 : arg_der_index.resize( arg_names.size() );
101 :
102 100 : if( argout.size()==0 ) {
103 48 : for(unsigned i=0; i<arg_names.size(); ++i) {
104 8 : argout.push_back( arg_names[i] );
105 8 : arg_der_index[i]=i;
106 : }
107 : } else {
108 88 : if(!disable_checks) {
109 88 : if( arg_names.size()!=argout.size() ) error("mismatched numbers of arguments in pdb frames");
110 : }
111 704 : for(unsigned i=0; i<arg_names.size(); ++i) {
112 176 : if(!disable_checks) {
113 176 : if( argout[i]!=arg_names[i] ) error("found mismatched arguments in pdb frames");
114 176 : arg_der_index[i]=i;
115 : } else {
116 : bool found=false;
117 0 : for(unsigned j=0; j<arg_names.size(); ++j) {
118 0 : if( argout[j]==arg_names[i] ) { found=true; arg_der_index[i]=j; break; }
119 : }
120 : if( !found ) {
121 0 : arg_der_index[i]=argout.size(); argout.push_back( arg_names[i] );
122 : }
123 : }
124 : }
125 : }
126 100 : }
127 :
128 265 : void ReferenceArguments::printArguments( OFile& ofile, const std::string& fmt ) const {
129 265 : if( arg_names.size()>0 ) {
130 263 : ofile.printf("REMARK ARG=%s", arg_names[0].c_str() );
131 1315 : for(unsigned i=1; i<arg_names.size(); ++i) ofile.printf(",%s", arg_names[i].c_str() );
132 263 : ofile.printf("\n");
133 :
134 263 : ofile.printf("REMARK ");
135 : std::string descr2;
136 263 : if(fmt.find("-")!=std::string::npos) {
137 0 : descr2="%s=" + fmt + " ";
138 : } else {
139 : // This ensures numbers are left justified (i.e. next to the equals sign
140 : std::size_t psign=fmt.find("%");
141 263 : plumed_assert( psign!=std::string::npos );
142 1052 : descr2="%s=%-" + fmt.substr(psign+1) + " ";
143 : }
144 2630 : for(unsigned i=0; i<arg_names.size(); ++i) ofile.printf( descr2.c_str(),arg_names[i].c_str(), reference_args[i] );
145 263 : ofile.printf("\n");
146 : }
147 : // Missing print out of metrics
148 265 : }
149 :
150 677 : const std::vector<double>& ReferenceArguments::getReferenceMetric() {
151 677 : if( hasmetric ) {
152 0 : unsigned ntot=(reference_args.size() / 2 )*(reference_args.size()+1);
153 0 : if( trig_metric.size()!=ntot ) trig_metric.resize( ntot );
154 : unsigned k=0;
155 0 : for(unsigned i=0; i<reference_args.size(); ++i) {
156 0 : for(unsigned j=i; j<reference_args.size(); ++j) {
157 : plumed_dbg_assert( fabs( metric(i,j)-metric(j,i) ) < epsilon );
158 0 : trig_metric[k]=metric(i,j); k++;
159 : }
160 : }
161 : } else {
162 677 : if( trig_metric.size()!=reference_args.size() ) trig_metric.resize( reference_args.size() );
163 5416 : for(unsigned i=0; i<reference_args.size(); ++i) trig_metric[i]=weights[i];
164 : }
165 677 : return trig_metric;
166 : }
167 :
168 48940 : double ReferenceArguments::calculateArgumentDistance( const std::vector<Value*> & vals, const std::vector<double>& arg,
169 : ReferenceValuePack& myder, const bool& squared ) const {
170 48940 : double r=0; std::vector<double> arg_ders( vals.size() );
171 48940 : if( hasmetric ) {
172 0 : for(unsigned i=0; i<reference_args.size(); ++i) {
173 0 : unsigned ik=arg_der_index[i]; arg_ders[ ik ]=0;
174 0 : double dp_i=vals[ik]->difference( reference_args[i], arg[ik] );
175 0 : for(unsigned j=0; j<reference_args.size(); ++j) {
176 : double dp_j;
177 0 : unsigned jk=arg_der_index[j];
178 0 : if(i==j) dp_j=dp_i;
179 0 : else dp_j=vals[jk]->difference( reference_args[j], arg[jk] );
180 :
181 0 : arg_ders[ ik ]+=2.0*metric(i,j)*dp_j; // Factor of two for off diagonal terms as you have terms from ij and ji
182 0 : r+=dp_i*dp_j*metric(i,j);
183 : }
184 : }
185 : } else {
186 391520 : for(unsigned i=0; i<reference_args.size(); ++i) {
187 97880 : unsigned ik=arg_der_index[i];
188 391520 : double dp_i=vals[ik]->difference( reference_args[i], arg[ik] );
189 195760 : r+=weights[i]*dp_i*dp_i; arg_ders[ik]=2.0*weights[i]*dp_i;
190 : }
191 : }
192 48940 : if(!squared) {
193 581 : r=sqrt(r); double ir=1.0/(2.0*r);
194 3486 : for(unsigned i=0; i<arg_ders.size(); ++i) myder.setArgumentDerivatives( i, arg_ders[i]*ir );
195 : } else {
196 290154 : for(unsigned i=0; i<arg_ders.size(); ++i) myder.setArgumentDerivatives( i, arg_ders[i] );
197 : }
198 48940 : return r;
199 : }
200 :
201 1964 : void ReferenceArguments::extractArgumentDisplacement( const std::vector<Value*>& vals, const std::vector<double>& arg, std::vector<double>& dirout ) const {
202 1964 : if( hasmetric ) {
203 0 : plumed_error();
204 : } else {
205 15712 : for(unsigned j=0; j<reference_args.size(); ++j) {
206 23568 : unsigned jk=arg_der_index[j]; dirout[jk]=sqrtweight[j]*vals[jk]->difference( reference_args[j], arg[jk] );
207 : }
208 : }
209 1964 : }
210 :
211 1294 : double ReferenceArguments::projectArgDisplacementOnVector( const std::vector<double>& eigv, const std::vector<Value*>& vals, const std::vector<double>& arg, ReferenceValuePack& mypack ) const {
212 1294 : if( hasmetric ) {
213 0 : plumed_error();
214 : } else {
215 : double proj=0;
216 10352 : for(unsigned j=0; j<reference_args.size(); ++j) {
217 2588 : unsigned jk=arg_der_index[j];
218 12940 : proj += eigv[j]*sqrtweight[j]*vals[jk]->difference( reference_args[j], arg[jk] );
219 5176 : mypack.setArgumentDerivatives( jk, eigv[j]*sqrtweight[j] );
220 : }
221 1294 : return proj;
222 : }
223 : }
224 :
225 500 : void ReferenceArguments::displaceReferenceArguments( const double& weight, const std::vector<double>& displace ) {
226 : plumed_dbg_assert( displace.size()==getNumberOfReferenceArguments() );
227 4864 : for(unsigned i=0; i<displace.size(); ++i) reference_args[i] += weight*displace[i];
228 500 : }
229 :
230 4839 : }
|
