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 "ReferenceConfiguration.h"
23 : #include "ReferenceArguments.h"
24 : #include "ReferenceAtoms.h"
25 : #include "Direction.h"
26 : #include "core/Value.h"
27 : #include "tools/OFile.h"
28 : #include "tools/PDB.h"
29 :
30 : namespace PLMD {
31 :
32 2207 : ReferenceConfigurationOptions::ReferenceConfigurationOptions( const std::string& type ):
33 2207 : tt(type)
34 : {
35 2207 : }
36 :
37 977 : bool ReferenceConfigurationOptions::usingFastOption() const {
38 1954 : return (tt.find("-FAST")!=std::string::npos);
39 : }
40 :
41 7 : std::string ReferenceConfigurationOptions::getMultiRMSDType() const {
42 14 : plumed_assert( tt.find("MULTI-")!=std::string::npos );
43 : std::size_t dot=tt.find_first_of("MULTI-");
44 7 : return tt.substr(dot+6);
45 : }
46 :
47 2207 : ReferenceConfiguration::ReferenceConfiguration( const ReferenceConfigurationOptions& ro ):
48 2207 : name(ro.tt)
49 : // arg_ders(0),
50 : // atom_ders(0)
51 : {
52 2207 : weight=0.0;
53 2207 : }
54 :
55 10684 : ReferenceConfiguration::~ReferenceConfiguration()
56 : {
57 2671 : }
58 :
59 723 : std::string ReferenceConfiguration::getName() const {
60 723 : return name;
61 : }
62 :
63 545 : void ReferenceConfiguration::set( const PDB& pdb ) {
64 545 : line=pdb.getRemark();
65 : std::string ignore;
66 1090 : if( parse("TYPE",ignore,true) ) {
67 48 : if(ignore!=name) error("mismatch for name");
68 : }
69 1090 : if( !parse("WEIGHT",weight,true) ) weight=1.0;
70 545 : read( pdb );
71 541 : }
72 :
73 : // void ReferenceConfiguration::setNumberOfArguments( const unsigned& n ){
74 : // arg_ders.resize(n); tmparg.resize(n);
75 : // }
76 :
77 : // void ReferenceConfiguration::setNumberOfAtoms( const unsigned& n ){
78 : // atom_ders.resize(n);
79 : // }
80 :
81 : // bool ReferenceConfiguration::getVirial( Tensor& virout ) const {
82 : // if(virialWasSet) virout=virial;
83 : // return virialWasSet;
84 : // }
85 :
86 18 : void ReferenceConfiguration::parseFlag( const std::string&key, bool&t ) {
87 18 : Tools::parseFlag(line,key,t);
88 18 : }
89 :
90 0 : void ReferenceConfiguration::error(const std::string& msg) {
91 0 : plumed_merror("error reading reference configuration of type " + name + " : " + msg );
92 : }
93 :
94 1614 : void ReferenceConfiguration::setNamesAndAtomNumbers( const std::vector<AtomNumber>& numbers, const std::vector<std::string>& arg ) {
95 1614 : ReferenceAtoms* atoms=dynamic_cast<ReferenceAtoms*>( this );
96 1614 : if(!atoms) {
97 857 : plumed_massert( numbers.size()==0, "expecting no atomic positions");
98 : //setNumberOfAtoms( 0 );
99 : } else {
100 757 : atoms->setAtomNumbers( numbers );
101 : // setNumberOfAtoms( numbers.size() );
102 : }
103 : // Copy the arguments to the reference
104 1614 : ReferenceArguments* args=dynamic_cast<ReferenceArguments*>( this );
105 1614 : if(!args) {
106 554 : plumed_massert( arg.size()==0, "expecting no arguments");
107 : // setNumberOfArguments(0);
108 : } else {
109 1060 : args->setArgumentNames( arg );
110 : // setNumberOfArguments( arg.size() );
111 : }
112 1614 : }
113 :
114 1988 : void ReferenceConfiguration::setReferenceConfig( const std::vector<Vector>& pos, const std::vector<double>& arg, const std::vector<double>& metric ) {
115 : // plumed_dbg_assert( pos.size()==atom_ders.size() && arg.size()==arg_ders.size() );
116 : // Copy the atomic positions to the reference
117 1988 : ReferenceAtoms* atoms=dynamic_cast<ReferenceAtoms*>( this );
118 1988 : if(!atoms) {
119 1424 : plumed_massert( pos.size()==0, "expecting no atomic positions");
120 : } else {
121 1692 : std::vector<double> align_in( pos.size(), 1.0 ), displace_in( pos.size(), 1.0 );
122 564 : atoms->setReferenceAtoms( pos, align_in, displace_in );
123 : }
124 : // Copy the arguments to the reference
125 1988 : ReferenceArguments* args=dynamic_cast<ReferenceArguments*>( this );
126 1988 : if(!args) {
127 1128 : plumed_massert( arg.size()==0 && metric.size()==0, "expecting no arguments");
128 : } else {
129 1424 : args->setReferenceArguments( arg, metric );
130 : }
131 1988 : }
132 :
133 448 : void ReferenceConfiguration::checkRead() {
134 448 : if(!line.empty()) {
135 0 : std::string msg="cannot understand the following words from the input line : ";
136 0 : for(unsigned i=0; i<line.size(); i++) {
137 0 : if(i>0) msg = msg + ", ";
138 0 : msg = msg + line[i];
139 : }
140 0 : error(msg);
141 : }
142 448 : }
143 :
144 174802 : double ReferenceConfiguration::calculate( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals,
145 : ReferenceValuePack& myder, const bool& squared ) const {
146 : // clearDerivatives();
147 174802 : std::vector<double> tmparg( vals.size() );
148 638196 : for(unsigned i=0; i<vals.size(); ++i) tmparg[i]=vals[i]->get();
149 349604 : return calc( pos, pbc, vals, tmparg, myder, squared );
150 : }
151 :
152 : // void ReferenceConfiguration::copyDerivatives( const ReferenceConfiguration* ref ){
153 : // plumed_dbg_assert( ref->atom_ders.size()==atom_ders.size() && ref->arg_ders.size()==arg_ders.size() );
154 : // for(unsigned i=0;i<atom_ders.size();++i) atom_ders[i]=ref->atom_ders[i];
155 : // for(unsigned i=0;i<arg_ders.size();++i) arg_ders[i]=ref->arg_ders[i];
156 : // virialWasSet=ref->virialWasSet; virial=ref->virial;
157 : // }
158 :
159 0 : void ReferenceConfiguration::print( OFile& ofile, const double& time, const double& weight, const double& lunits, const double& old_norm ) {
160 0 : ofile.printf("REMARK TIME=%f LOG_WEIGHT=%f OLD_NORM=%f\n",time, weight, old_norm );
161 0 : print( ofile, "%f", lunits ); // HARD CODED FORMAT HERE AS THIS IS FOR CHECKPOINT FILE
162 0 : }
163 :
164 279 : void ReferenceConfiguration::print( OFile& ofile, const std::string& fmt, const double& lunits ) {
165 837 : ofile.printf("REMARK TYPE=%s\n",getName().c_str() );
166 279 : ReferenceArguments* args=dynamic_cast<ReferenceArguments*>(this);
167 279 : if(args) args->printArguments( ofile, fmt );
168 279 : ReferenceAtoms* atoms=dynamic_cast<ReferenceAtoms*>(this);
169 279 : if(atoms) atoms->printAtoms( ofile, lunits );
170 279 : ofile.printf("END\n");
171 279 : }
172 :
173 500 : void ReferenceConfiguration::displaceReferenceConfiguration( const double& weight, Direction& dir ) {
174 500 : ReferenceArguments* args=dynamic_cast<ReferenceArguments*>(this);
175 500 : if( args ) args->displaceReferenceArguments( weight, dir.getReferenceArguments() );
176 500 : ReferenceAtoms* atoms=dynamic_cast<ReferenceAtoms*>(this);
177 500 : if( atoms ) atoms->displaceReferenceAtoms( weight, dir.getReferencePositions() );
178 500 : }
179 :
180 3079 : void ReferenceConfiguration::extractDisplacementVector( const std::vector<Vector>& pos, const std::vector<Value*>& vals,
181 : const std::vector<double>& arg, const bool& nflag,
182 : Direction& mydir ) const {
183 3079 : const ReferenceAtoms* atoms=dynamic_cast<const ReferenceAtoms*>( this );
184 3079 : if( atoms ) atoms->extractAtomicDisplacement( pos, mydir.reference_atoms );
185 3079 : const ReferenceArguments* args=dynamic_cast<const ReferenceArguments*>( this );
186 3079 : if( args ) args->extractArgumentDisplacement( vals, arg, mydir.reference_args );
187 :
188 : // Normalize direction if required
189 3079 : if( nflag ) {
190 : // Calculate length of vector
191 8 : double tmp, norm=0; mydir.normalized = true;
192 184 : for(unsigned i=0; i<mydir.getReferencePositions().size(); ++i) {
193 392 : for(unsigned k=0; k<3; ++k) { tmp=mydir.getReferencePositions()[i][k]; norm+=tmp*tmp; }
194 : }
195 16 : for(unsigned i=0; i<mydir.getReferenceArguments().size(); ++i) { tmp=mydir.getReferenceArguments()[i]; norm+=tmp*tmp; }
196 8 : norm = sqrt( norm );
197 : // And normalize
198 184 : for(unsigned i=0; i<mydir.getReferencePositions().size(); ++i) {
199 392 : for(unsigned k=0; k<3; ++k) { mydir.reference_atoms[i][k] /=norm; }
200 : }
201 16 : for(unsigned i=0; i<mydir.getReferenceArguments().size(); ++i) { mydir.reference_args[i] /= norm; }
202 : }
203 3079 : }
204 :
205 2474 : double ReferenceConfiguration::projectDisplacementOnVector( const Direction& mydir,
206 : const std::vector<Value*>& vals, const std::vector<double>& arg,
207 : ReferenceValuePack& mypack ) const {
208 : double proj=0;
209 2474 : const ReferenceAtoms* atoms=dynamic_cast<const ReferenceAtoms*>( this );
210 2474 : if( atoms ) proj += atoms->projectAtomicDisplacementOnVector( mydir.normalized, mydir.getReferencePositions(), mypack );
211 2474 : const ReferenceArguments* args=dynamic_cast<const ReferenceArguments*>( this );
212 2474 : if( args ) proj += args->projectArgDisplacementOnVector( mydir.getReferenceArguments(), vals, arg, mypack );
213 2474 : return proj;
214 : }
215 :
216 9900 : double distance( const Pbc& pbc, const std::vector<Value*> & vals, ReferenceConfiguration* ref1, ReferenceConfiguration* ref2, const bool& squared ) {
217 : unsigned nder;
218 19800 : if( ref1->getReferencePositions().size()>0 ) nder=ref1->getReferenceArguments().size() + 3*ref1->getReferencePositions().size() + 9;
219 19800 : else nder=ref1->getReferenceArguments().size();
220 :
221 39600 : MultiValue myvals( 1, nder ); ReferenceValuePack myder( ref1->getReferenceArguments().size(), ref1->getReferencePositions().size(), myvals );
222 9900 : double dist1=ref1->calc( ref2->getReferencePositions(), pbc, vals, ref2->getReferenceArguments(), myder, squared );
223 : #ifndef NDEBUG
224 : // Check that A - B = B - A
225 : double dist2=ref2->calc( ref1->getReferencePositions(), pbc, vals, ref1->getReferenceArguments(), myder, squared );
226 : plumed_dbg_assert( fabs(dist1-dist2)<epsilon );
227 : #endif
228 9900 : return dist1;
229 : }
230 :
231 4839 : }
|
