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 "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             : #include "core/GenericMolInfo.h"
      30             : 
      31             : namespace PLMD {
      32             : 
      33      519115 : ReferenceConfigurationOptions::ReferenceConfigurationOptions( const std::string& type ):
      34      519115 :   tt(type) {
      35      519115 : }
      36             : 
      37         455 : bool ReferenceConfigurationOptions::usingFastOption() const {
      38         455 :   return (tt.find("-FAST")!=std::string::npos);
      39             : }
      40             : 
      41           5 : std::string ReferenceConfigurationOptions::getMultiRMSDType() const {
      42           5 :   plumed_assert( tt.find("MULTI-")!=std::string::npos );
      43           5 :   std::size_t dot=tt.find_first_of("MULTI-");
      44           5 :   return tt.substr(dot+6);
      45             : }
      46             : 
      47      519115 : ReferenceConfiguration::ReferenceConfiguration( const ReferenceConfigurationOptions& ro ):
      48      519115 :   name(ro.tt) {
      49      519115 : }
      50             : 
      51      519596 : ReferenceConfiguration::~ReferenceConfiguration() {
      52     1039192 : }
      53             : 
      54         284 : std::string ReferenceConfiguration::getName() const {
      55         284 :   return name;
      56             : }
      57             : 
      58           0 : [[noreturn]] void ReferenceConfiguration::error(const std::string& msg) {
      59           0 :   plumed_merror("error reading reference configuration of type " + name + " : " + msg );
      60             : }
      61             : 
      62      175359 : double ReferenceConfiguration::calculate( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals,
      63             :     ReferenceValuePack& myder, const bool& squared ) const {
      64      175359 :   std::vector<double> tmparg( vals.size() );
      65      248599 :   for(unsigned i=0; i<vals.size(); ++i) {
      66       73240 :     tmparg[i]=vals[i]->get();
      67             :   }
      68      350718 :   return calc( pos, pbc, vals, tmparg, myder, squared );
      69             : }
      70             : 
      71         500 : void ReferenceConfiguration::displaceReferenceConfiguration( const double& weight, Direction& dir ) {
      72         500 :   ReferenceArguments* args=dynamic_cast<ReferenceArguments*>(this);
      73         500 :   if( args ) {
      74         500 :     args->displaceReferenceArguments( weight, dir.getReferenceArguments() );
      75             :   }
      76         500 :   ReferenceAtoms* atoms=dynamic_cast<ReferenceAtoms*>(this);
      77         500 :   if( atoms ) {
      78         494 :     atoms->displaceReferenceAtoms( weight, dir.getReferencePositions() );
      79             :   }
      80         500 : }
      81             : 
      82        3081 : void ReferenceConfiguration::extractDisplacementVector( const std::vector<Vector>& pos, const std::vector<Value*>& vals,
      83             :     const std::vector<double>& arg, const bool& nflag,
      84             :     Direction& mydir ) const {
      85        3081 :   const ReferenceAtoms* atoms=dynamic_cast<const ReferenceAtoms*>( this );
      86        3081 :   if( atoms ) {
      87        1117 :     atoms->extractAtomicDisplacement( pos, mydir.reference_atoms );
      88             :   }
      89        3081 :   const ReferenceArguments* args=dynamic_cast<const ReferenceArguments*>( this );
      90        3081 :   if( args ) {
      91        1974 :     args->extractArgumentDisplacement( vals, arg, mydir.reference_args );
      92             :   }
      93             : 
      94             :   // Normalize direction if required
      95        3081 :   if( nflag ) {
      96             :     // Calculate length of vector
      97             :     double tmp, norm=0;
      98          10 :     mydir.normalized = true;
      99          80 :     for(unsigned i=0; i<mydir.getReferencePositions().size(); ++i) {
     100         280 :       for(unsigned k=0; k<3; ++k) {
     101         210 :         tmp=mydir.getReferencePositions()[i][k];
     102         210 :         norm+=tmp*tmp;
     103             :       }
     104             :     }
     105          10 :     for(unsigned i=0; i<mydir.getReferenceArguments().size(); ++i) {
     106           0 :       tmp=mydir.getReferenceArguments()[i];
     107           0 :       norm+=tmp*tmp;
     108             :     }
     109          10 :     norm = std::sqrt( norm );
     110             :     // And normalize
     111          80 :     for(unsigned i=0; i<mydir.getReferencePositions().size(); ++i) {
     112         280 :       for(unsigned k=0; k<3; ++k) {
     113         210 :         mydir.reference_atoms[i][k] /=norm;
     114             :       }
     115             :     }
     116          10 :     for(unsigned i=0; i<mydir.getReferenceArguments().size(); ++i) {
     117           0 :       mydir.reference_args[i] /= norm;
     118             :     }
     119             :   }
     120        3081 : }
     121             : 
     122        3588 : double ReferenceConfiguration::projectDisplacementOnVector( const Direction& mydir,
     123             :     const std::vector<Value*>& vals, const std::vector<double>& arg,
     124             :     ReferenceValuePack& mypack ) const {
     125             :   double proj=0;
     126        3588 :   const ReferenceAtoms* atoms=dynamic_cast<const ReferenceAtoms*>( this );
     127        3588 :   if( atoms ) {
     128        1202 :     proj += atoms->projectAtomicDisplacementOnVector( mydir.normalized, mydir.getReferencePositions(), mypack );
     129             :   }
     130        3588 :   const ReferenceArguments* args=dynamic_cast<const ReferenceArguments*>( this );
     131        3588 :   if( args ) {
     132        2386 :     proj += args->projectArgDisplacementOnVector( mydir.getReferenceArguments(), vals, arg, mypack );
     133             :   }
     134        3588 :   return proj;
     135             : }
     136             : 
     137      259099 : double distance( const Pbc& pbc, const std::vector<Value*> & vals, ReferenceConfiguration* ref1, ReferenceConfiguration* ref2, const bool& squared ) {
     138             :   unsigned nder;
     139      259099 :   if( ref1->getReferencePositions().size()>0 ) {
     140           7 :     nder=ref1->getReferenceArguments().size() + 3*ref1->getReferencePositions().size() + 9;
     141             :   } else {
     142      259092 :     nder=ref1->getReferenceArguments().size();
     143             :   }
     144             : 
     145      259099 :   MultiValue myvals( 1, nder );
     146      259099 :   ReferenceValuePack myder( ref1->getReferenceArguments().size(), ref1->getReferencePositions().size(), myvals );
     147      259099 :   double dist1=ref1->calc( ref2->getReferencePositions(), pbc, vals, ref2->getReferenceArguments(), myder, squared );
     148             : #ifndef NDEBUG
     149             :   // Check that A - B = B - A
     150             :   double dist2=ref2->calc( ref1->getReferencePositions(), pbc, vals, ref1->getReferenceArguments(), myder, squared );
     151             :   plumed_dbg_assert( std::fabs(dist1-dist2)<epsilon );
     152             : #endif
     153      259099 :   return dist1;
     154      259099 : }
     155             : 
     156             : 
     157             : }