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 13 : void Mapping::registerKeywords( Keywords& keys ) {
35 13 : Action::registerKeywords( keys );
36 13 : ActionWithValue::registerKeywords( keys );
37 13 : ActionWithArguments::registerKeywords( keys );
38 13 : ActionAtomistic::registerKeywords( keys );
39 13 : vesselbase::ActionWithVessel::registerKeywords( keys );
40 26 : keys.add("compulsory","REFERENCE","a pdb file containing the set of reference configurations");
41 26 : keys.add("compulsory","PROPERTY","the property to be used in the index. This should be in the REMARK of the reference");
42 26 : 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 26 : keys.addFlag("DISABLE_CHECKS",false,"disable checks on reference input structures.");
46 13 : }
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 : {
55 : // Read the input
56 10 : std::string mtype; parse("TYPE",mtype);
57 10 : bool skipchecks; parseFlag("DISABLE_CHECKS",skipchecks);
58 :
59 : // Read the properties we require
60 : bool ispath=false;
61 20 : if( keywords.exists("PROPERTY") ) {
62 6 : std::vector<std::string> propnames; parseVector("PROPERTY",propnames);
63 3 : if(propnames.size()==0) error("no properties were specified");
64 9 : for(unsigned i=0; i<propnames.size(); ++i) property.insert( std::pair<std::string,std::vector<double> >( propnames[i], std::vector<double>() ) );
65 3 : } else {
66 14 : property.insert( std::pair<std::string,std::vector<double> >( "spath", std::vector<double>() ) ); ispath=true;
67 : }
68 :
69 : // Open reference file
70 10 : std::string reference; parse("REFERENCE",reference);
71 10 : FILE* fp=fopen(reference.c_str(),"r");
72 10 : if(!fp) error("could not open reference file " + reference );
73 :
74 : // Read all reference configurations
75 : bool do_read=true; unsigned nfram=0; double wnorm=0., ww;
76 426 : while (do_read) {
77 : // Read the pdb file
78 852 : PDB mypdb; do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
79 : // Break if we are done
80 426 : if( !do_read ) break ;
81 : // Check for required properties
82 416 : if( !ispath ) {
83 : double prop;
84 378 : for(std::map<std::string,std::vector<double> >::iterator it=property.begin(); it!=property.end(); ++it) {
85 252 : if( !mypdb.getArgumentValue( it->first, prop ) ) error("pdb input does not have contain property named " + it->first );
86 252 : it->second.push_back(prop);
87 : }
88 : } else {
89 580 : property.find("spath")->second.push_back( myframes.size()+1 );
90 : }
91 : // Fix argument names
92 416 : expandArgKeywordInPDB( mypdb );
93 : // And read the frame
94 416 : myframes.emplace_back( metricRegister().create<ReferenceConfiguration>( mtype, mypdb ) );
95 832 : if( !mypdb.getArgumentValue( "WEIGHT", ww ) ) ww=1.0;
96 416 : weights.push_back( ww ); wnorm+=ww; nfram++;
97 426 : }
98 10 : fclose(fp);
99 :
100 10 : if(nfram==0 ) error("no reference configurations were specified");
101 10 : log.printf(" found %u configurations in file %s\n",nfram,reference.c_str() );
102 426 : for(unsigned i=0; i<weights.size(); ++i) weights[i] = weights[i]/wnorm;
103 :
104 : // Finish the setup of the mapping object
105 : // Get the arguments and atoms that are required
106 : std::vector<AtomNumber> atoms; std::vector<std::string> args;
107 426 : for(unsigned i=0; i<myframes.size(); ++i) { myframes[i]->getAtomRequests( atoms, skipchecks ); myframes[i]->getArgumentRequests( args, skipchecks ); }
108 10 : std::vector<Value*> req_args; interpretArgumentList( args, req_args );
109 10 : if( req_args.size()>0 && atoms.size()>0 ) error("cannot mix atoms and arguments");
110 10 : if( req_args.size()>0 ) requestArguments( req_args );
111 10 : if( atoms.size()>0 ) {
112 8 : log.printf(" found %zu atoms in input \n",atoms.size());
113 8 : log.printf(" with indices : ");
114 112 : for(unsigned i=0; i<atoms.size(); ++i) {
115 104 : if(i%25==0) log<<"\n";
116 104 : log.printf("%d ",atoms[i].serial());
117 : }
118 8 : log.printf("\n");
119 8 : requestAtoms( atoms );
120 : }
121 : // Duplicate all frames (duplicates are used by sketch-map)
122 : // mymap->duplicateFrameList();
123 : // fframes.resize( 2*nfram, 0.0 ); dfframes.resize( 2*nfram, 0.0 );
124 : // plumed_assert( !mymap->mappingNeedsSetup() );
125 : // Resize all derivative arrays
126 : // mymap->setNumberOfAtomsAndArguments( atoms.size(), args.size() );
127 : // Resize forces array
128 10 : if( getNumberOfAtoms()>0 ) {
129 8 : forcesToApply.resize( 3*getNumberOfAtoms() + 9 + getNumberOfArguments() );
130 : } else {
131 2 : forcesToApply.resize( getNumberOfArguments() );
132 : }
133 20 : }
134 :
135 19 : void Mapping::turnOnDerivatives() {
136 19 : ActionWithValue::turnOnDerivatives();
137 19 : needsDerivatives();
138 19 : }
139 :
140 0 : double Mapping::getLambda() {
141 0 : plumed_merror("lambda is not defined in this mapping type");
142 : }
143 :
144 0 : std::string Mapping::getArgumentName( unsigned& iarg ) {
145 0 : if( iarg < getNumberOfArguments() ) return getPntrToArgument(iarg)->getName();
146 0 : unsigned iatom=iarg - getNumberOfArguments();
147 0 : std::string atnum; Tools::convert( getAbsoluteIndex(iatom).serial(),atnum);
148 0 : unsigned icomp=iatom%3;
149 0 : if(icomp==0) return "pos" + atnum + "x";
150 0 : if(icomp==1) return "pos" + atnum + "y";
151 0 : return "pos" + atnum + "z";
152 : }
153 :
154 172372 : void Mapping::finishPackSetup( const unsigned& ifunc, ReferenceValuePack& mypack ) const {
155 : mypack.setValIndex(0);
156 172372 : unsigned nargs2=myframes[ifunc]->getNumberOfReferenceArguments();
157 172372 : unsigned nat2=myframes[ifunc]->getNumberOfReferencePositions();
158 172372 : if( mypack.getNumberOfAtoms()!=nat2 || mypack.getNumberOfArguments()!=nargs2 ) mypack.resize( nargs2, nat2 );
159 172372 : if( nat2>0 ) {
160 137592 : ReferenceAtoms* myat2=dynamic_cast<ReferenceAtoms*>( myframes[ifunc].get() ); plumed_dbg_assert( myat2 );
161 1926288 : for(unsigned i=0; i<nat2; ++i) mypack.setAtomIndex( i, myat2->getAtomIndex(i) );
162 : }
163 172372 : }
164 :
165 172372 : double Mapping::calculateDistanceFunction( const unsigned& ifunc, ReferenceValuePack& myder, const bool& squared ) const {
166 : // Calculate the distance
167 172372 : double dd = myframes[ifunc]->calculate( getPositions(), getPbc(), getArguments(), myder, squared );
168 : // Transform distance by whatever
169 172372 : double df, ff=transformHD( dd, df ); myder.scaleAllDerivatives( df );
170 : // And the virial
171 172372 : if( getNumberOfAtoms()>0 && !myder.virialWasSet() ) {
172 137592 : Tensor tvir; tvir.zero();
173 1926288 : for(unsigned i=0; i<myder.getNumberOfAtoms(); ++i) tvir +=-1.0*Tensor( getPosition( myder.getAtomIndex(i) ), myder.getAtomDerivative(i) );
174 137592 : myder.addBoxDerivatives( tvir );
175 : }
176 172372 : return ff;
177 : }
178 :
179 13175 : ReferenceConfiguration* Mapping::getReferenceConfiguration( const unsigned& ifunc ) {
180 13175 : return myframes[ifunc].get();
181 : }
182 :
183 0 : void Mapping::calculateNumericalDerivatives( ActionWithValue* a ) {
184 0 : if( getNumberOfArguments()>0 ) {
185 0 : ActionWithArguments::calculateNumericalDerivatives( a );
186 : }
187 0 : if( getNumberOfAtoms()>0 ) {
188 0 : Matrix<double> save_derivatives( getNumberOfComponents(), getNumberOfArguments() );
189 0 : for(int j=0; j<getNumberOfComponents(); ++j) {
190 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) save_derivatives(j,i)=getPntrToComponent(j)->getDerivative(i);
191 : }
192 0 : calculateAtomicNumericalDerivatives( a, getNumberOfArguments() );
193 0 : for(int j=0; j<getNumberOfComponents(); ++j) {
194 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) getPntrToComponent(j)->addDerivative( i, save_derivatives(j,i) );
195 : }
196 : }
197 0 : }
198 :
199 5015 : void Mapping::apply() {
200 5015 : if( getForcesFromVessels( forcesToApply ) ) {
201 0 : addForcesOnArguments( forcesToApply );
202 0 : if( getNumberOfAtoms()>0 ) setForcesOnAtoms( forcesToApply, getNumberOfArguments() );
203 : }
204 5015 : }
205 :
206 : }
207 : }
208 :
209 :
|
