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