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 "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 9 : void Mapping::registerKeywords( Keywords& keys ) {
35 9 : Action::registerKeywords( keys );
36 9 : ActionWithValue::registerKeywords( keys );
37 9 : ActionWithArguments::registerKeywords( keys );
38 9 : ActionAtomistic::registerKeywords( keys );
39 9 : vesselbase::ActionWithVessel::registerKeywords( keys );
40 9 : keys.add("compulsory","REFERENCE","a pdb file containing the set of reference configurations");
41 9 : keys.add("compulsory","PROPERTY","the property to be used in the index. This should be in the REMARK of the reference");
42 : 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 9 : "\\ref dists");
45 9 : keys.addFlag("DISABLE_CHECKS",false,"disable checks on reference input structures.");
46 9 : }
47 :
48 7 : Mapping::Mapping(const ActionOptions&ao):
49 : Action(ao),
50 : ActionAtomistic(ao),
51 : ActionWithArguments(ao),
52 : ActionWithValue(ao),
53 7 : ActionWithVessel(ao)
54 : {
55 : // Read the input
56 7 : std::string mtype; parse("TYPE",mtype);
57 7 : bool skipchecks; parseFlag("DISABLE_CHECKS",skipchecks);
58 : // Setup the object that does the mapping
59 7 : mymap = new PointWiseMapping( mtype, skipchecks );
60 :
61 : // Read the properties we require
62 7 : if( keywords.exists("PROPERTY") ) {
63 3 : std::vector<std::string> property;
64 3 : parseVector("PROPERTY",property);
65 3 : if(property.size()==0) error("no properties were specified");
66 3 : mymap->setPropertyNames( property, false );
67 : } else {
68 4 : std::vector<std::string> property(1);
69 4 : property[0]="spath";
70 4 : mymap->setPropertyNames( property, true );
71 : }
72 :
73 : // Open reference file
74 14 : std::string reference; parse("REFERENCE",reference);
75 7 : FILE* fp=fopen(reference.c_str(),"r");
76 7 : if(!fp) error("could not open reference file " + reference );
77 :
78 : // Read all reference configurations
79 14 : bool do_read=true; std::vector<double> weights;
80 7 : unsigned nfram=0, wnorm=0., ww;
81 308 : while (do_read) {
82 301 : PDB mypdb;
83 : // Read the pdb file
84 301 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
85 : // Fix argument names
86 301 : expandArgKeywordInPDB( mypdb );
87 301 : if(do_read) {
88 294 : mymap->readFrame( mypdb ); ww=mymap->getWeight( nfram );
89 294 : weights.push_back( ww );
90 294 : wnorm+=ww; nfram++;
91 : } else {
92 7 : break;
93 : }
94 294 : }
95 7 : fclose(fp);
96 :
97 7 : if(nfram==0 ) error("no reference configurations were specified");
98 7 : log.printf(" found %u configurations in file %s\n",nfram,reference.c_str() );
99 7 : for(unsigned i=0; i<weights.size(); ++i) weights[i] /= wnorm;
100 7 : mymap->setWeights( weights );
101 :
102 : // Finish the setup of the mapping object
103 : // Get the arguments and atoms that are required
104 14 : std::vector<AtomNumber> atoms; std::vector<std::string> args;
105 7 : mymap->getAtomAndArgumentRequirements( atoms, args );
106 14 : std::vector<Value*> req_args; interpretArgumentList( args, req_args );
107 7 : if( req_args.size()>0 && atoms.size()>0 ) error("cannot mix atoms and arguments");
108 7 : if( req_args.size()>0 ) requestArguments( req_args );
109 7 : if( atoms.size()>0 ) requestAtoms( atoms );
110 : // Duplicate all frames (duplicates are used by sketch-map)
111 7 : mymap->duplicateFrameList();
112 : // fframes.resize( 2*nfram, 0.0 ); dfframes.resize( 2*nfram, 0.0 );
113 7 : plumed_assert( !mymap->mappingNeedsSetup() );
114 : // Resize all derivative arrays
115 : // mymap->setNumberOfAtomsAndArguments( atoms.size(), args.size() );
116 : // Resize forces array
117 7 : if( getNumberOfAtoms()>0 ) {
118 7 : forcesToApply.resize( 3*getNumberOfAtoms() + 9 + getNumberOfArguments() );
119 : } else {
120 0 : forcesToApply.resize( getNumberOfArguments() );
121 7 : }
122 7 : }
123 :
124 15 : void Mapping::turnOnDerivatives() {
125 15 : ActionWithValue::turnOnDerivatives();
126 15 : needsDerivatives();
127 15 : }
128 :
129 14 : Mapping::~Mapping() {
130 7 : delete mymap;
131 7 : }
132 :
133 14 : unsigned Mapping::getPropertyIndex( const std::string& name ) const {
134 14 : return mymap->getPropertyIndex( name );
135 : }
136 :
137 0 : double Mapping::getLambda() {
138 0 : plumed_merror("lambda is not defined in this mapping type");
139 : }
140 :
141 0 : std::string Mapping::getArgumentName( unsigned& iarg ) {
142 0 : if( iarg < getNumberOfArguments() ) return getPntrToArgument(iarg)->getName();
143 0 : unsigned iatom=iarg - getNumberOfArguments();
144 0 : std::string atnum; Tools::convert( getAbsoluteIndex(iatom).serial(),atnum);
145 0 : unsigned icomp=iatom%3;
146 0 : if(icomp==0) return "pos" + atnum + "x";
147 0 : if(icomp==1) return "pos" + atnum + "y";
148 0 : return "pos" + atnum + "z";
149 : }
150 :
151 112808 : void Mapping::finishPackSetup( const unsigned& ifunc, ReferenceValuePack& mypack ) const {
152 112808 : ReferenceConfiguration* myref=mymap->getFrame(ifunc); mypack.setValIndex(0);
153 113481 : unsigned nargs2=myref->getNumberOfReferenceArguments(); unsigned nat2=myref->getNumberOfReferencePositions();
154 114428 : if( mypack.getNumberOfAtoms()!=nat2 || mypack.getNumberOfArguments()!=nargs2 ) mypack.resize( nargs2, nat2 );
155 114447 : if( nat2>0 ) {
156 113484 : ReferenceAtoms* myat2=dynamic_cast<ReferenceAtoms*>( myref ); plumed_dbg_assert( myat2 );
157 113717 : for(unsigned i=0; i<nat2; ++i) mypack.setAtomIndex( i, myat2->getAtomIndex(i) );
158 : }
159 115538 : }
160 :
161 113910 : double Mapping::calculateDistanceFunction( const unsigned& ifunc, ReferenceValuePack& myder, const bool& squared ) const {
162 : // Calculate the distance
163 113910 : double dd = mymap->calcDistanceFromConfiguration( ifunc, getPositions(), getPbc(), getArguments(), myder, squared );
164 : // Transform distance by whatever
165 114122 : double df, ff=transformHD( dd, df ); myder.scaleAllDerivatives( df );
166 : // And the virial
167 114534 : if( getNumberOfAtoms()>0 && !myder.virialWasSet() ) {
168 114410 : Tensor tvir; tvir.zero();
169 114201 : for(unsigned i=0; i<myder.getNumberOfAtoms(); ++i) tvir +=-1.0*Tensor( getPosition( myder.getAtomIndex(i) ), myder.getAtomDerivative(i) );
170 114658 : myder.addBoxDerivatives( tvir );
171 : }
172 114649 : return ff;
173 : }
174 :
175 0 : ReferenceConfiguration* Mapping::getReferenceConfiguration( const unsigned& ifunc ) {
176 0 : return mymap->getFrame( ifunc );
177 : }
178 :
179 0 : void Mapping::calculateNumericalDerivatives( ActionWithValue* a ) {
180 0 : if( getNumberOfArguments()>0 ) {
181 0 : ActionWithArguments::calculateNumericalDerivatives( a );
182 : }
183 0 : if( getNumberOfAtoms()>0 ) {
184 0 : Matrix<double> save_derivatives( getNumberOfComponents(), getNumberOfArguments() );
185 0 : for(unsigned j=0; j<getNumberOfComponents(); ++j) {
186 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) save_derivatives(j,i)=getPntrToComponent(j)->getDerivative(i);
187 : }
188 0 : calculateAtomicNumericalDerivatives( a, getNumberOfArguments() );
189 0 : for(unsigned j=0; j<getNumberOfComponents(); ++j) {
190 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) getPntrToComponent(j)->addDerivative( i, save_derivatives(j,i) );
191 0 : }
192 : }
193 0 : }
194 :
195 3822 : void Mapping::apply() {
196 3822 : if( getForcesFromVessels( forcesToApply ) ) {
197 0 : addForcesOnArguments( forcesToApply );
198 0 : if( getNumberOfAtoms()>0 ) setForcesOnAtoms( forcesToApply, getNumberOfArguments() );
199 : }
200 3822 : }
201 :
202 : }
203 2523 : }
204 :
205 :
|
