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 "ReferenceArguments.h"
      23             : #include "ReferenceAtoms.h"
      24             : #include "tools/OFile.h"
      25             : #include "core/Value.h"
      26             : #include "tools/PDB.h"
      27             : 
      28             : namespace PLMD {
      29             : 
      30      518583 : ReferenceArguments::ReferenceArguments( const ReferenceConfigurationOptions& ro ):
      31             :   ReferenceConfiguration(ro),
      32      518583 :   hasweights(false),
      33      518583 :   hasmetric(false)
      34             : {
      35      518583 : }
      36             : 
      37      519050 : void ReferenceArguments::readArgumentsFromPDB( const PDB& pdb ) {
      38      519050 :   ReferenceAtoms* aref=dynamic_cast<ReferenceAtoms*>( this );
      39      519050 :   arg_names.resize( pdb.getArgumentNames().size() );
      40     2015570 :   for(unsigned i=0; i<arg_names.size(); ++i) arg_names[i]=pdb.getArgumentNames()[i];
      41      519050 :   if( !aref && arg_names.size()==0 ) error("no arguments in input PDB file");
      42             : 
      43      519050 :   reference_args.resize( arg_names.size() ); arg_der_index.resize( arg_names.size() );
      44     2015570 :   for(unsigned i=0; i<arg_names.size(); ++i) {
      45     1496520 :     if( !pdb.getArgumentValue(arg_names[i], reference_args[i]) ) error("argument " + arg_names[i] + " was not set in pdb input");
      46     1496520 :     arg_der_index[i]=i;
      47             :   }
      48             : 
      49      519050 :   if( hasweights ) {
      50           0 :     plumed_massert( !hasmetric, "should not have weights if we are using metric");
      51           0 :     weights.resize( arg_names.size() ); sqrtweight.resize( arg_names.size() );
      52           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
      53           0 :       if( !pdb.getArgumentValue("sigma_" + arg_names[i], weights[i]) ) error("value sigma_" + arg_names[i] + " was not set in pdb input");
      54           0 :       sqrtweight[i] = std::sqrt( weights[i] );
      55             :     }
      56      519050 :   } else if( hasmetric ) {
      57             :     plumed_massert( !hasweights, "should not have weights if we are using metric");
      58           0 :     double thissig; metric.resize( arg_names.size(), arg_names.size() );
      59           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
      60           0 :       for(unsigned j=i; j<reference_args.size(); ++j) {
      61           0 :         if( !pdb.getArgumentValue("sigma_" + arg_names[i] + "_" + arg_names[j], thissig) ) {
      62           0 :           error("value sigma_" + arg_names[i] + "_" + arg_names[j] + " was not set in pdb input");
      63             :         }
      64           0 :         metric(i,j)=metric(j,i)=thissig;
      65             :       }
      66             :     }
      67             :   } else {
      68      519050 :     weights.resize( arg_names.size() ); sqrtweight.resize( arg_names.size() );
      69     2015570 :     for(unsigned i=0; i<weights.size(); ++i) sqrtweight[i]=weights[i]=1.0;
      70             :   }
      71      519050 : }
      72             : 
      73         494 : void ReferenceArguments::setReferenceArguments( const std::vector<double>& arg_vals, const std::vector<double>& sigma ) {
      74         494 :   moveReferenceArguments( arg_vals );
      75             : 
      76         494 :   if( hasmetric ) {
      77             :     unsigned k=0;
      78           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
      79           0 :       for(unsigned j=i; j<reference_args.size(); ++j) {
      80           0 :         metric(i,j)=metric(j,i)=sigma[k]; k++;
      81             :       }
      82             :     }
      83           0 :     plumed_assert( k==sigma.size() );
      84             :   } else {
      85         494 :     plumed_assert( reference_args.size()==sigma.size() );
      86        1448 :     for(unsigned i=0; i<reference_args.size(); ++i) weights[i]=sigma[i];
      87             :   }
      88         494 : }
      89             : 
      90         494 : void ReferenceArguments::moveReferenceArguments( const std::vector<double>& arg_vals ) {
      91             :   plumed_dbg_assert( reference_args.size()==arg_vals.size() );
      92        1448 :   for(unsigned i=0; i<arg_vals.size(); ++i) reference_args[i]=arg_vals[i];
      93         494 : }
      94             : 
      95         187 : void ReferenceArguments::getArgumentRequests( std::vector<std::string>& argout, bool disable_checks ) {
      96         187 :   arg_der_index.resize( arg_names.size() );
      97             : 
      98         187 :   if( argout.size()==0 ) {
      99          31 :     for(unsigned i=0; i<arg_names.size(); ++i) {
     100          15 :       argout.push_back( arg_names[i] );
     101          15 :       arg_der_index[i]=i;
     102             :     }
     103             :   } else {
     104         171 :     if(!disable_checks) {
     105         171 :       if( arg_names.size()!=argout.size() ) error("mismatched numbers of arguments in pdb frames");
     106             :     }
     107         762 :     for(unsigned i=0; i<arg_names.size(); ++i) {
     108         591 :       if(!disable_checks) {
     109         591 :         if( argout[i]!=arg_names[i] ) error("found mismatched arguments in pdb frames");
     110         591 :         arg_der_index[i]=i;
     111             :       } else {
     112             :         bool found=false;
     113           0 :         for(unsigned j=0; j<arg_names.size(); ++j) {
     114           0 :           if( argout[j]==arg_names[i] ) { found=true; arg_der_index[i]=j; break; }
     115             :         }
     116             :         if( !found ) {
     117           0 :           arg_der_index[i]=argout.size(); argout.push_back( arg_names[i] );
     118             :         }
     119             :       }
     120             :     }
     121             :   }
     122         187 : }
     123             : 
     124           0 : const std::vector<double>& ReferenceArguments::getReferenceMetric() {
     125           0 :   if( hasmetric ) {
     126           0 :     unsigned ntot=(reference_args.size() / 2 )*(reference_args.size()+1);
     127           0 :     if( trig_metric.size()!=ntot ) trig_metric.resize( ntot );
     128             :     unsigned k=0;
     129           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
     130           0 :       for(unsigned j=i; j<reference_args.size(); ++j) {
     131             :         plumed_dbg_assert( std::fabs( metric(i,j)-metric(j,i) ) < epsilon );
     132           0 :         trig_metric[k]=metric(i,j); k++;
     133             :       }
     134             :     }
     135             :   } else {
     136           0 :     if( trig_metric.size()!=reference_args.size() ) trig_metric.resize( reference_args.size() );
     137           0 :     for(unsigned i=0; i<reference_args.size(); ++i) trig_metric[i]=weights[i];
     138             :   }
     139           0 :   return trig_metric;
     140             : }
     141             : 
     142      299760 : double ReferenceArguments::calculateArgumentDistance( const std::vector<Value*> & vals, const std::vector<double>& arg,
     143             :     ReferenceValuePack& myder, const bool& squared ) const {
     144      299760 :   double r=0; std::vector<double> arg_ders( vals.size() );
     145      299760 :   if( hasmetric ) {
     146           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
     147           0 :       unsigned ik=arg_der_index[i]; arg_ders[ ik ]=0;
     148           0 :       double dp_i=vals[ik]->difference( reference_args[i], arg[ik] );
     149           0 :       for(unsigned j=0; j<reference_args.size(); ++j) {
     150             :         double dp_j;
     151           0 :         unsigned jk=arg_der_index[j];
     152           0 :         if(i==j) dp_j=dp_i;
     153           0 :         else dp_j=vals[jk]->difference( reference_args[j], arg[jk] );
     154             : 
     155           0 :         arg_ders[ ik ]+=2.0*metric(i,j)*dp_j;    // Factor of two for off diagonal terms as you have terms from ij and ji
     156           0 :         r+=dp_i*dp_j*metric(i,j);
     157             :       }
     158             :     }
     159             :   } else {
     160     1128420 :     for(unsigned i=0; i<reference_args.size(); ++i) {
     161      828660 :       unsigned ik=arg_der_index[i];
     162      828660 :       double dp_i=vals[ik]->difference( reference_args[i], arg[ik] );
     163      828660 :       r+=weights[i]*dp_i*dp_i; arg_ders[ik]=2.0*weights[i]*dp_i;
     164             :     }
     165             :   }
     166      299760 :   if(!squared) {
     167         572 :     r=std::sqrt(r); double ir=1.0/(2.0*r);
     168        1716 :     for(unsigned i=0; i<arg_ders.size(); ++i) myder.setArgumentDerivatives( i, arg_ders[i]*ir );
     169             :   } else {
     170     1126704 :     for(unsigned i=0; i<arg_ders.size(); ++i) myder.setArgumentDerivatives( i, arg_ders[i] );
     171             :   }
     172      299760 :   return r;
     173             : }
     174             : 
     175        1964 : void ReferenceArguments::extractArgumentDisplacement( const std::vector<Value*>& vals, const std::vector<double>& arg, std::vector<double>& dirout ) const {
     176        1964 :   if( hasmetric ) {
     177           0 :     plumed_error();
     178             :   } else {
     179        5892 :     for(unsigned j=0; j<reference_args.size(); ++j) {
     180        3928 :       unsigned jk=arg_der_index[j]; dirout[jk]=sqrtweight[j]*vals[jk]->difference( reference_args[j], arg[jk] );
     181             :     }
     182             :   }
     183        1964 : }
     184             : 
     185        2386 : double ReferenceArguments::projectArgDisplacementOnVector( const std::vector<double>& eigv, const std::vector<Value*>& vals, const std::vector<double>& arg, ReferenceValuePack& mypack ) const {
     186        2386 :   if( hasmetric ) {
     187           0 :     plumed_error();
     188             :   } else {
     189             :     double proj=0;
     190        7158 :     for(unsigned j=0; j<reference_args.size(); ++j) {
     191        4772 :       unsigned jk=arg_der_index[j];
     192        4772 :       proj += eigv[j]*sqrtweight[j]*vals[jk]->difference( reference_args[j], arg[jk] );
     193        4772 :       mypack.setArgumentDerivatives( jk, eigv[j]*sqrtweight[j] );
     194             :     }
     195        2386 :     return proj;
     196             :   }
     197             : }
     198             : 
     199         500 : void ReferenceArguments::displaceReferenceArguments( const double& weight, const std::vector<double>& displace ) {
     200             :   plumed_dbg_assert( displace.size()==getNumberOfReferenceArguments() );
     201        1466 :   for(unsigned i=0; i<displace.size(); ++i) reference_args[i] += weight*displace[i];
     202         500 : }
     203             : 
     204             : }