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 "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 52 : keys.add("compulsory","REFERENCE","a pdb file containing the set of reference configurations");
41 52 : keys.add("compulsory","PROPERTY","the property to be used in the index. This should be in the REMARK of the reference");
42 65 : 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 39 : 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 20 : ActionWithVessel(ao)
54 : {
55 : // Read the input
56 20 : std::string mtype; parse("TYPE",mtype);
57 20 : bool skipchecks; parseFlag("DISABLE_CHECKS",skipchecks);
58 : // Setup the object that does the mapping
59 10 : mymap = new PointWiseMapping( mtype, skipchecks );
60 :
61 : // Read the properties we require
62 20 : if( keywords.exists("PROPERTY") ) {
63 3 : std::vector<std::string> property;
64 6 : parseVector("PROPERTY",property);
65 3 : if(property.size()==0) error("no properties were specified");
66 3 : mymap->setPropertyNames( property, false );
67 : } else {
68 14 : std::vector<std::string> property(1);
69 : property[0]="spath";
70 7 : mymap->setPropertyNames( property, true );
71 : }
72 :
73 : // Open reference file
74 20 : std::string reference; parse("REFERENCE",reference);
75 20 : FILE* fp=fopen(reference.c_str(),"r");
76 10 : if(!fp) error("could not open reference file " + reference );
77 :
78 : // Read all reference configurations
79 : bool do_read=true; std::vector<double> weights;
80 : unsigned nfram=0, wnorm=0., ww;
81 426 : while (do_read) {
82 842 : PDB mypdb;
83 : // Read the pdb file
84 852 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
85 : // Fix argument names
86 426 : expandArgKeywordInPDB( mypdb );
87 426 : if(do_read) {
88 832 : mymap->readFrame( mypdb ); ww=mymap->getWeight( nfram );
89 832 : weights.push_back( ww );
90 416 : wnorm+=ww; nfram++;
91 : } else {
92 : break;
93 : }
94 : }
95 10 : fclose(fp);
96 :
97 10 : if(nfram==0 ) error("no reference configurations were specified");
98 20 : log.printf(" found %u configurations in file %s\n",nfram,reference.c_str() );
99 1268 : for(unsigned i=0; i<weights.size(); ++i) weights[i] /= wnorm;
100 10 : mymap->setWeights( weights );
101 :
102 : // Finish the setup of the mapping object
103 : // Get the arguments and atoms that are required
104 10 : std::vector<AtomNumber> atoms; std::vector<std::string> args;
105 10 : mymap->getAtomAndArgumentRequirements( atoms, args );
106 10 : std::vector<Value*> req_args; interpretArgumentList( args, req_args );
107 12 : if( req_args.size()>0 && atoms.size()>0 ) error("cannot mix atoms and arguments");
108 10 : if( req_args.size()>0 ) requestArguments( req_args );
109 10 : if( atoms.size()>0 ) requestAtoms( atoms );
110 : // Duplicate all frames (duplicates are used by sketch-map)
111 : // mymap->duplicateFrameList();
112 : // fframes.resize( 2*nfram, 0.0 ); dfframes.resize( 2*nfram, 0.0 );
113 : // plumed_assert( !mymap->mappingNeedsSetup() );
114 : // Resize all derivative arrays
115 : // mymap->setNumberOfAtomsAndArguments( atoms.size(), args.size() );
116 : // Resize forces array
117 : if( getNumberOfAtoms()>0 ) {
118 16 : forcesToApply.resize( 3*getNumberOfAtoms() + 9 + getNumberOfArguments() );
119 : } else {
120 2 : forcesToApply.resize( getNumberOfArguments() );
121 : }
122 10 : }
123 :
124 19 : void Mapping::turnOnDerivatives() {
125 19 : ActionWithValue::turnOnDerivatives();
126 19 : needsDerivatives();
127 19 : }
128 :
129 30 : Mapping::~Mapping() {
130 10 : delete mymap;
131 10 : }
132 :
133 14 : unsigned Mapping::getPropertyIndex( const std::string& name ) const {
134 14 : return mymap->getPropertyIndex( name );
135 : }
136 :
137 20 : void Mapping::setPropertyValue( const unsigned& iframe, const unsigned& jprop, const double& property ) {
138 20 : mymap->setProjectionCoordinate( iframe, jprop, property );
139 20 : }
140 :
141 0 : double Mapping::getLambda() {
142 0 : plumed_merror("lambda is not defined in this mapping type");
143 : }
144 :
145 0 : std::string Mapping::getArgumentName( unsigned& iarg ) {
146 0 : if( iarg < getNumberOfArguments() ) return getPntrToArgument(iarg)->getName();
147 0 : unsigned iatom=iarg - getNumberOfArguments();
148 0 : std::string atnum; Tools::convert( getAbsoluteIndex(iatom).serial(),atnum);
149 0 : unsigned icomp=iatom%3;
150 0 : if(icomp==0) return "pos" + atnum + "x";
151 0 : if(icomp==1) return "pos" + atnum + "y";
152 0 : return "pos" + atnum + "z";
153 : }
154 :
155 172372 : void Mapping::finishPackSetup( const unsigned& ifunc, ReferenceValuePack& mypack ) const {
156 172372 : ReferenceConfiguration* myref=mymap->getFrame(ifunc); mypack.setValIndex(0);
157 172372 : unsigned nargs2=myref->getNumberOfReferenceArguments(); unsigned nat2=myref->getNumberOfReferencePositions();
158 344744 : if( mypack.getNumberOfAtoms()!=nat2 || mypack.getNumberOfArguments()!=nargs2 ) mypack.resize( nargs2, nat2 );
159 172372 : if( nat2>0 ) {
160 137592 : ReferenceAtoms* myat2=dynamic_cast<ReferenceAtoms*>( myref ); plumed_dbg_assert( myat2 );
161 3714984 : 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 344744 : double dd = mymap->calcDistanceFromConfiguration( ifunc, 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 309964 : if( getNumberOfAtoms()>0 && !myder.virialWasSet() ) {
172 137592 : Tensor tvir; tvir.zero();
173 5641272 : 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 26350 : return mymap->getFrame( ifunc );
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(unsigned 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(unsigned 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 4839 : }
208 :
209 :
|
