Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2023 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 "Mapping.h"
23 : #include "vesselbase/Vessel.h"
24 : #include "reference/MetricRegister.h"
25 : #include "reference/ReferenceAtoms.h"
26 : #include "tools/PDB.h"
27 : #include "tools/Matrix.h"
28 : #include "core/PlumedMain.h"
29 : #include "core/Atoms.h"
30 :
31 : namespace PLMD {
32 : namespace mapping {
33 :
34 22 : void Mapping::registerKeywords( Keywords& keys ) {
35 22 : Action::registerKeywords( keys );
36 22 : ActionWithValue::registerKeywords( keys );
37 22 : ActionWithArguments::registerKeywords( keys );
38 22 : ActionAtomistic::registerKeywords( keys );
39 22 : vesselbase::ActionWithVessel::registerKeywords( keys );
40 44 : keys.add("compulsory","REFERENCE","a pdb file containing the set of reference configurations");
41 44 : keys.add("compulsory","PROPERTY","the property to be used in the index. This should be in the REMARK of the reference");
42 44 : keys.add("compulsory","TYPE","OPTIMAL-FAST","the manner in which distances are calculated. More information on the different "
43 : "metrics that are available in PLUMED can be found in the section of the manual on "
44 : "\\ref dists");
45 44 : keys.addFlag("DISABLE_CHECKS",false,"disable checks on reference input structures.");
46 22 : }
47 :
48 10 : Mapping::Mapping(const ActionOptions&ao):
49 : Action(ao),
50 : ActionAtomistic(ao),
51 : ActionWithArguments(ao),
52 : ActionWithValue(ao),
53 10 : ActionWithVessel(ao) {
54 : // Read the input
55 : std::string mtype;
56 10 : parse("TYPE",mtype);
57 : bool skipchecks;
58 10 : parseFlag("DISABLE_CHECKS",skipchecks);
59 :
60 : // Read the properties we require
61 : bool ispath=false;
62 20 : if( keywords.exists("PROPERTY") ) {
63 : std::vector<std::string> propnames;
64 6 : parseVector("PROPERTY",propnames);
65 3 : if(propnames.size()==0) {
66 0 : error("no properties were specified");
67 : }
68 9 : for(unsigned i=0; i<propnames.size(); ++i) {
69 12 : property.insert( std::pair<std::string,std::vector<double> >( propnames[i], std::vector<double>() ) );
70 : }
71 3 : } else {
72 14 : property.insert( std::pair<std::string,std::vector<double> >( "spath", std::vector<double>() ) );
73 : ispath=true;
74 : }
75 :
76 : // Open reference file
77 : std::string reference;
78 10 : parse("REFERENCE",reference);
79 :
80 10 : FILE* fp=this->fopen(reference.c_str(),"r");
81 10 : if(!fp) {
82 0 : error("could not open reference file " + reference );
83 : }
84 :
85 : // call fclose when exiting this block
86 10 : auto deleter=[this](FILE* f) {
87 10 : this->fclose(f);
88 10 : };
89 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
90 :
91 : // Read all reference configurations
92 : bool do_read=true;
93 : unsigned nfram=0;
94 : double wnorm=0., ww;
95 426 : while (do_read) {
96 : // Read the pdb file
97 426 : PDB mypdb;
98 852 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
99 : // Break if we are done
100 426 : if( !do_read ) {
101 : break ;
102 : }
103 : // Check for required properties
104 416 : if( !ispath ) {
105 : double prop;
106 378 : for(std::map<std::string,std::vector<double> >::iterator it=property.begin(); it!=property.end(); ++it) {
107 252 : if( !mypdb.getArgumentValue( it->first, prop ) ) {
108 0 : error("pdb input does not have contain property named " + it->first );
109 : }
110 252 : it->second.push_back(prop);
111 : }
112 : } else {
113 580 : property.find("spath")->second.push_back( myframes.size()+1 );
114 : }
115 : // Fix argument names
116 416 : expandArgKeywordInPDB( mypdb );
117 : // And read the frame
118 416 : myframes.emplace_back( metricRegister().create<ReferenceConfiguration>( mtype, mypdb ) );
119 832 : if( !mypdb.getArgumentValue( "WEIGHT", ww ) ) {
120 416 : ww=1.0;
121 : }
122 416 : weights.push_back( ww );
123 416 : wnorm+=ww;
124 416 : nfram++;
125 426 : }
126 :
127 10 : if(nfram==0 ) {
128 0 : error("no reference configurations were specified");
129 : }
130 10 : log.printf(" found %u configurations in file %s\n",nfram,reference.c_str() );
131 426 : for(unsigned i=0; i<weights.size(); ++i) {
132 416 : weights[i] = weights[i]/wnorm;
133 : }
134 :
135 : // Finish the setup of the mapping object
136 : // Get the arguments and atoms that are required
137 : std::vector<AtomNumber> atoms;
138 : std::vector<std::string> args;
139 426 : for(unsigned i=0; i<myframes.size(); ++i) {
140 416 : myframes[i]->getAtomRequests( atoms, skipchecks );
141 416 : myframes[i]->getArgumentRequests( args, skipchecks );
142 : }
143 : std::vector<Value*> req_args;
144 10 : interpretArgumentList( args, req_args );
145 10 : if( req_args.size()>0 && atoms.size()>0 ) {
146 0 : error("cannot mix atoms and arguments");
147 : }
148 10 : if( req_args.size()>0 ) {
149 2 : requestArguments( req_args );
150 : }
151 10 : if( atoms.size()>0 ) {
152 8 : log.printf(" found %zu atoms in input \n",atoms.size());
153 8 : log.printf(" with indices : ");
154 112 : for(unsigned i=0; i<atoms.size(); ++i) {
155 104 : if(i%25==0) {
156 8 : log<<"\n";
157 : }
158 104 : log.printf("%d ",atoms[i].serial());
159 : }
160 8 : log.printf("\n");
161 8 : requestAtoms( atoms );
162 : }
163 : // Duplicate all frames (duplicates are used by sketch-map)
164 : // mymap->duplicateFrameList();
165 : // fframes.resize( 2*nfram, 0.0 ); dfframes.resize( 2*nfram, 0.0 );
166 : // plumed_assert( !mymap->mappingNeedsSetup() );
167 : // Resize all derivative arrays
168 : // mymap->setNumberOfAtomsAndArguments( atoms.size(), args.size() );
169 : // Resize forces array
170 10 : if( getNumberOfAtoms()>0 ) {
171 8 : forcesToApply.resize( 3*getNumberOfAtoms() + 9 + getNumberOfArguments() );
172 : } else {
173 2 : forcesToApply.resize( getNumberOfArguments() );
174 : }
175 30 : }
176 :
177 19 : void Mapping::turnOnDerivatives() {
178 19 : ActionWithValue::turnOnDerivatives();
179 19 : needsDerivatives();
180 19 : }
181 :
182 0 : double Mapping::getLambda() {
183 0 : plumed_merror("lambda is not defined in this mapping type");
184 : }
185 :
186 0 : std::string Mapping::getArgumentName( unsigned& iarg ) {
187 0 : if( iarg < getNumberOfArguments() ) {
188 0 : return getPntrToArgument(iarg)->getName();
189 : }
190 0 : unsigned iatom=iarg - getNumberOfArguments();
191 : std::string atnum;
192 0 : Tools::convert( getAbsoluteIndex(iatom).serial(),atnum);
193 0 : unsigned icomp=iatom%3;
194 0 : if(icomp==0) {
195 0 : return "pos" + atnum + "x";
196 : }
197 0 : if(icomp==1) {
198 0 : return "pos" + atnum + "y";
199 : }
200 0 : return "pos" + atnum + "z";
201 : }
202 :
203 172372 : void Mapping::finishPackSetup( const unsigned& ifunc, ReferenceValuePack& mypack ) const {
204 : mypack.setValIndex(0);
205 172372 : unsigned nargs2=myframes[ifunc]->getNumberOfReferenceArguments();
206 172372 : unsigned nat2=myframes[ifunc]->getNumberOfReferencePositions();
207 172372 : if( mypack.getNumberOfAtoms()!=nat2 || mypack.getNumberOfArguments()!=nargs2 ) {
208 0 : mypack.resize( nargs2, nat2 );
209 : }
210 172372 : if( nat2>0 ) {
211 137592 : ReferenceAtoms* myat2=dynamic_cast<ReferenceAtoms*>( myframes[ifunc].get() );
212 : plumed_dbg_assert( myat2 );
213 1926288 : for(unsigned i=0; i<nat2; ++i) {
214 : mypack.setAtomIndex( i, myat2->getAtomIndex(i) );
215 : }
216 : }
217 172372 : }
218 :
219 172372 : double Mapping::calculateDistanceFunction( const unsigned& ifunc, ReferenceValuePack& myder, const bool& squared ) const {
220 : // Calculate the distance
221 172372 : double dd = myframes[ifunc]->calculate( getPositions(), getPbc(), getArguments(), myder, squared );
222 : // Transform distance by whatever
223 172372 : double df, ff=transformHD( dd, df );
224 172372 : myder.scaleAllDerivatives( df );
225 : // And the virial
226 172372 : if( getNumberOfAtoms()>0 && !myder.virialWasSet() ) {
227 137592 : Tensor tvir;
228 137592 : tvir.zero();
229 1926288 : for(unsigned i=0; i<myder.getNumberOfAtoms(); ++i) {
230 1788696 : tvir +=-1.0*Tensor( getPosition( myder.getAtomIndex(i) ), myder.getAtomDerivative(i) );
231 : }
232 137592 : myder.addBoxDerivatives( tvir );
233 : }
234 172372 : return ff;
235 : }
236 :
237 13175 : ReferenceConfiguration* Mapping::getReferenceConfiguration( const unsigned& ifunc ) {
238 13175 : return myframes[ifunc].get();
239 : }
240 :
241 0 : void Mapping::calculateNumericalDerivatives( ActionWithValue* a ) {
242 0 : if( getNumberOfArguments()>0 ) {
243 0 : ActionWithArguments::calculateNumericalDerivatives( a );
244 : }
245 0 : if( getNumberOfAtoms()>0 ) {
246 0 : Matrix<double> save_derivatives( getNumberOfComponents(), getNumberOfArguments() );
247 0 : for(int j=0; j<getNumberOfComponents(); ++j) {
248 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
249 0 : save_derivatives(j,i)=getPntrToComponent(j)->getDerivative(i);
250 : }
251 : }
252 0 : calculateAtomicNumericalDerivatives( a, getNumberOfArguments() );
253 0 : for(int j=0; j<getNumberOfComponents(); ++j) {
254 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
255 0 : getPntrToComponent(j)->addDerivative( i, save_derivatives(j,i) );
256 : }
257 : }
258 : }
259 0 : }
260 :
261 5015 : void Mapping::apply() {
262 5015 : if( getForcesFromVessels( forcesToApply ) ) {
263 0 : addForcesOnArguments( forcesToApply );
264 0 : if( getNumberOfAtoms()>0 ) {
265 0 : setForcesOnAtoms( forcesToApply, getNumberOfArguments() );
266 : }
267 : }
268 5015 : }
269 :
270 : }
271 : }
272 :
273 :
|
