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 "ReferenceArguments.h"
23 : #include "ReferenceAtoms.h"
24 : #include "tools/OFile.h"
25 : #include "core/Value.h"
26 :
27 : namespace PLMD {
28 :
29 859 : ReferenceArguments::ReferenceArguments( const ReferenceConfigurationOptions& ro ):
30 : ReferenceConfiguration(ro),
31 : hasweights(false),
32 859 : hasmetric(false)
33 : {
34 859 : }
35 :
36 8 : void ReferenceArguments::readArgumentsFromPDB( const PDB& pdb ) {
37 8 : ReferenceAtoms* aref=dynamic_cast<ReferenceAtoms*>( this );
38 8 : if( !aref ) parseVector( "ARG", arg_names );
39 8 : else parseVector( "ARG", arg_names, true );
40 :
41 8 : reference_args.resize( arg_names.size() );
42 8 : for(unsigned i=0; i<arg_names.size(); ++i) parse( arg_names[i], reference_args[i] );
43 :
44 8 : if( hasweights ) {
45 0 : plumed_massert( !hasmetric, "should not have weights if we are using metric");
46 0 : weights.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 : }
50 8 : } else if( hasmetric ) {
51 0 : plumed_massert( !hasweights, "should not have weights if we are using metric");
52 0 : double thissig; metric.resize( arg_names.size(), arg_names.size() );
53 0 : for(unsigned i=0; i<reference_args.size(); ++i) {
54 0 : for(unsigned j=i; j<reference_args.size(); ++j) {
55 0 : parse( "sigma_" + arg_names[i] + "_" + arg_names[j], thissig );
56 0 : metric(i,j)=metric(j,i)=thissig;
57 : }
58 : }
59 : } else {
60 8 : weights.resize( arg_names.size() );
61 8 : for(unsigned i=0; i<weights.size(); ++i) weights[i]=1.0;
62 : }
63 8 : }
64 :
65 849 : void ReferenceArguments::setArgumentNames( const std::vector<std::string>& arg_vals ) {
66 849 : reference_args.resize( arg_vals.size() );
67 849 : arg_names.resize( arg_vals.size() );
68 849 : arg_der_index.resize( arg_vals.size() );
69 2545 : for(unsigned i=0; i<arg_vals.size(); ++i) {
70 1696 : arg_names[i]=arg_vals[i]; arg_der_index[i]=i;
71 : }
72 849 : if( hasmetric ) metric.resize( arg_vals.size(), arg_vals.size() );
73 849 : else weights.resize( arg_vals.size() );
74 849 : }
75 :
76 951 : void ReferenceArguments::setReferenceArguments( const std::vector<double>& arg_vals, const std::vector<double>& sigma ) {
77 : plumed_dbg_assert( reference_args.size()==arg_vals.size() );
78 951 : for(unsigned i=0; i<arg_vals.size(); ++i) reference_args[i]=arg_vals[i];
79 :
80 951 : if( hasmetric ) {
81 0 : 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 951 : plumed_assert( reference_args.size()==sigma.size() );
90 951 : for(unsigned i=0; i<reference_args.size(); ++i) weights[i]=sigma[i];
91 : }
92 951 : }
93 :
94 8 : void ReferenceArguments::getArgumentRequests( std::vector<std::string>& argout, bool disable_checks ) {
95 8 : arg_der_index.resize( arg_names.size() );
96 :
97 8 : if( argout.size()==0 ) {
98 8 : for(unsigned i=0; i<arg_names.size(); ++i) {
99 0 : argout.push_back( arg_names[i] );
100 0 : arg_der_index[i]=i;
101 : }
102 : } else {
103 0 : if(!disable_checks) {
104 0 : if( arg_names.size()!=argout.size() ) error("mismatched numbers of arguments in pdb frames");
105 : }
106 0 : for(unsigned i=0; i<arg_names.size(); ++i) {
107 0 : bool found=false;
108 0 : if(!disable_checks) {
109 0 : if( argout[i]!=arg_names[i] ) error("found mismatched arguments in pdb frames");
110 0 : arg_der_index[i]=i;
111 : } else {
112 0 : for(unsigned j=0; j<arg_names.size(); ++j) {
113 0 : if( argout[j]==arg_names[i] ) { found=true; arg_der_index[i]=j; break; }
114 : }
115 0 : if( !found ) {
116 0 : arg_der_index[i]=argout.size(); argout.push_back( arg_names[i] );
117 : }
118 : }
119 : }
120 : }
121 8 : }
122 :
123 205 : void ReferenceArguments::printArguments( OFile& ofile, const std::string& fmt ) const {
124 205 : if( arg_names.size()>0 ) {
125 203 : ofile.printf("REMARK ARG=%s", arg_names[0].c_str() );
126 203 : for(unsigned i=1; i<arg_names.size(); ++i) ofile.printf(",%s", arg_names[i].c_str() );
127 203 : ofile.printf("\n");
128 :
129 203 : ofile.printf("REMARK ");
130 203 : std::string descr2;
131 203 : if(fmt.find("-")!=std::string::npos) {
132 0 : descr2="%s=" + fmt + " ";
133 : } else {
134 : // This ensures numbers are left justified (i.e. next to the equals sign
135 203 : std::size_t psign=fmt.find("%");
136 203 : plumed_assert( psign!=std::string::npos );
137 203 : descr2="%s=%-" + fmt.substr(psign+1) + " ";
138 : }
139 203 : for(unsigned i=0; i<arg_names.size(); ++i) ofile.printf( descr2.c_str(),arg_names[i].c_str(), reference_args[i] );
140 203 : ofile.printf("\n");
141 : }
142 : // Missing print out of metrics
143 205 : }
144 :
145 200 : const std::vector<double>& ReferenceArguments::getReferenceMetric() {
146 200 : if( hasmetric ) {
147 0 : unsigned ntot=(reference_args.size() / 2 )*(reference_args.size()+1);
148 0 : if( trig_metric.size()!=ntot ) trig_metric.resize( ntot );
149 0 : unsigned k=0;
150 0 : for(unsigned i=0; i<reference_args.size(); ++i) {
151 0 : for(unsigned j=i; j<reference_args.size(); ++j) {
152 : plumed_dbg_assert( fabs( metric(i,j)-metric(j,i) ) < epsilon );
153 0 : trig_metric[k]=metric(i,j); k++;
154 : }
155 : }
156 : } else {
157 200 : if( trig_metric.size()!=reference_args.size() ) trig_metric.resize( reference_args.size() );
158 200 : for(unsigned i=0; i<reference_args.size(); ++i) trig_metric[i]=weights[i];
159 : }
160 200 : return trig_metric;
161 : }
162 :
163 10991 : double ReferenceArguments::calculateArgumentDistance( const std::vector<Value*> & vals, const std::vector<double>& arg,
164 : ReferenceValuePack& myder, const bool& squared ) const {
165 10991 : double r=0; std::vector<double> arg_ders( vals.size() );
166 10991 : if( hasmetric ) {
167 0 : for(unsigned i=0; i<reference_args.size(); ++i) {
168 0 : unsigned ik=arg_der_index[i]; arg_ders[ ik ]=0;
169 0 : double dp_i=vals[ik]->difference( reference_args[i], arg[ik] );
170 0 : for(unsigned j=0; j<reference_args.size(); ++j) {
171 : double dp_j;
172 0 : unsigned jk=arg_der_index[j];
173 0 : if(i==j) dp_j=dp_i;
174 0 : else dp_j=vals[jk]->difference( reference_args[j], arg[jk] );
175 :
176 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
177 0 : r+=dp_i*dp_j*metric(i,j);
178 : }
179 : }
180 : } else {
181 32973 : for(unsigned i=0; i<reference_args.size(); ++i) {
182 21982 : unsigned ik=arg_der_index[i];
183 21982 : double dp_i=vals[ik]->difference( reference_args[i], arg[ik] );
184 21982 : r+=weights[i]*dp_i*dp_i; arg_ders[ik]=2.0*weights[i]*dp_i;
185 : }
186 : }
187 10991 : if(!squared) {
188 0 : r=sqrt(r); double ir=1.0/(2.0*r);
189 0 : for(unsigned i=0; i<arg_ders.size(); ++i) myder.setArgumentDerivatives( i, arg_ders[i]*ir );
190 : } else {
191 10991 : for(unsigned i=0; i<arg_ders.size(); ++i) myder.setArgumentDerivatives( i, arg_ders[i] );
192 : }
193 10991 : return r;
194 : }
195 2523 : }
|
