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      518583 : }
      35             : 
      36      519050 : void ReferenceArguments::readArgumentsFromPDB( const PDB& pdb ) {
      37      519050 :   ReferenceAtoms* aref=dynamic_cast<ReferenceAtoms*>( this );
      38      519050 :   arg_names.resize( pdb.getArgumentNames().size() );
      39     2015570 :   for(unsigned i=0; i<arg_names.size(); ++i) {
      40     2993040 :     arg_names[i]=pdb.getArgumentNames()[i];
      41             :   }
      42      519050 :   if( !aref && arg_names.size()==0 ) {
      43           0 :     error("no arguments in input PDB file");
      44             :   }
      45             : 
      46      519050 :   reference_args.resize( arg_names.size() );
      47      519050 :   arg_der_index.resize( arg_names.size() );
      48     2015570 :   for(unsigned i=0; i<arg_names.size(); ++i) {
      49     1496520 :     if( !pdb.getArgumentValue(arg_names[i], reference_args[i]) ) {
      50           0 :       error("argument " + arg_names[i] + " was not set in pdb input");
      51             :     }
      52     1496520 :     arg_der_index[i]=i;
      53             :   }
      54             : 
      55      519050 :   if( hasweights ) {
      56           0 :     plumed_massert( !hasmetric, "should not have weights if we are using metric");
      57           0 :     weights.resize( arg_names.size() );
      58           0 :     sqrtweight.resize( arg_names.size() );
      59           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
      60           0 :       if( !pdb.getArgumentValue("sigma_" + arg_names[i], weights[i]) ) {
      61           0 :         error("value sigma_" + arg_names[i] + " was not set in pdb input");
      62             :       }
      63           0 :       sqrtweight[i] = std::sqrt( weights[i] );
      64             :     }
      65      519050 :   } else if( hasmetric ) {
      66             :     plumed_massert( !hasweights, "should not have weights if we are using metric");
      67             :     double thissig;
      68           0 :     metric.resize( arg_names.size(), arg_names.size() );
      69           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
      70           0 :       for(unsigned j=i; j<reference_args.size(); ++j) {
      71           0 :         if( !pdb.getArgumentValue("sigma_" + arg_names[i] + "_" + arg_names[j], thissig) ) {
      72           0 :           error("value sigma_" + arg_names[i] + "_" + arg_names[j] + " was not set in pdb input");
      73             :         }
      74           0 :         metric(i,j)=metric(j,i)=thissig;
      75             :       }
      76             :     }
      77             :   } else {
      78      519050 :     weights.resize( arg_names.size() );
      79      519050 :     sqrtweight.resize( arg_names.size() );
      80     2015570 :     for(unsigned i=0; i<weights.size(); ++i) {
      81     1496520 :       sqrtweight[i]=weights[i]=1.0;
      82             :     }
      83             :   }
      84      519050 : }
      85             : 
      86         494 : void ReferenceArguments::setReferenceArguments( const std::vector<double>& arg_vals, const std::vector<double>& sigma ) {
      87         494 :   moveReferenceArguments( arg_vals );
      88             : 
      89         494 :   if( hasmetric ) {
      90             :     unsigned k=0;
      91           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
      92           0 :       for(unsigned j=i; j<reference_args.size(); ++j) {
      93           0 :         metric(i,j)=metric(j,i)=sigma[k];
      94           0 :         k++;
      95             :       }
      96             :     }
      97           0 :     plumed_assert( k==sigma.size() );
      98             :   } else {
      99         494 :     plumed_assert( reference_args.size()==sigma.size() );
     100        1448 :     for(unsigned i=0; i<reference_args.size(); ++i) {
     101         954 :       weights[i]=sigma[i];
     102             :     }
     103             :   }
     104         494 : }
     105             : 
     106         494 : void ReferenceArguments::moveReferenceArguments( const std::vector<double>& arg_vals ) {
     107             :   plumed_dbg_assert( reference_args.size()==arg_vals.size() );
     108        1448 :   for(unsigned i=0; i<arg_vals.size(); ++i) {
     109         954 :     reference_args[i]=arg_vals[i];
     110             :   }
     111         494 : }
     112             : 
     113         187 : void ReferenceArguments::getArgumentRequests( std::vector<std::string>& argout, bool disable_checks ) {
     114         187 :   arg_der_index.resize( arg_names.size() );
     115             : 
     116         187 :   if( argout.size()==0 ) {
     117          31 :     for(unsigned i=0; i<arg_names.size(); ++i) {
     118          15 :       argout.push_back( arg_names[i] );
     119          15 :       arg_der_index[i]=i;
     120             :     }
     121             :   } else {
     122         171 :     if(!disable_checks) {
     123         171 :       if( arg_names.size()!=argout.size() ) {
     124           0 :         error("mismatched numbers of arguments in pdb frames");
     125             :       }
     126             :     }
     127         762 :     for(unsigned i=0; i<arg_names.size(); ++i) {
     128         591 :       if(!disable_checks) {
     129         591 :         if( argout[i]!=arg_names[i] ) {
     130           0 :           error("found mismatched arguments in pdb frames");
     131             :         }
     132         591 :         arg_der_index[i]=i;
     133             :       } else {
     134             :         bool found=false;
     135           0 :         for(unsigned j=0; j<arg_names.size(); ++j) {
     136           0 :           if( argout[j]==arg_names[i] ) {
     137             :             found=true;
     138           0 :             arg_der_index[i]=j;
     139             :             break;
     140             :           }
     141             :         }
     142             :         if( !found ) {
     143           0 :           arg_der_index[i]=argout.size();
     144           0 :           argout.push_back( arg_names[i] );
     145             :         }
     146             :       }
     147             :     }
     148             :   }
     149         187 : }
     150             : 
     151           0 : const std::vector<double>& ReferenceArguments::getReferenceMetric() {
     152           0 :   if( hasmetric ) {
     153           0 :     unsigned ntot=(reference_args.size() / 2 )*(reference_args.size()+1);
     154           0 :     if( trig_metric.size()!=ntot ) {
     155           0 :       trig_metric.resize( ntot );
     156             :     }
     157             :     unsigned k=0;
     158           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
     159           0 :       for(unsigned j=i; j<reference_args.size(); ++j) {
     160             :         plumed_dbg_assert( std::fabs( metric(i,j)-metric(j,i) ) < epsilon );
     161           0 :         trig_metric[k]=metric(i,j);
     162           0 :         k++;
     163             :       }
     164             :     }
     165             :   } else {
     166           0 :     if( trig_metric.size()!=reference_args.size() ) {
     167           0 :       trig_metric.resize( reference_args.size() );
     168             :     }
     169           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
     170           0 :       trig_metric[i]=weights[i];
     171             :     }
     172             :   }
     173           0 :   return trig_metric;
     174             : }
     175             : 
     176      299760 : double ReferenceArguments::calculateArgumentDistance( const std::vector<Value*> & vals, const std::vector<double>& arg,
     177             :     ReferenceValuePack& myder, const bool& squared ) const {
     178             :   double r=0;
     179      299760 :   std::vector<double> arg_ders( vals.size() );
     180      299760 :   if( hasmetric ) {
     181           0 :     for(unsigned i=0; i<reference_args.size(); ++i) {
     182           0 :       unsigned ik=arg_der_index[i];
     183           0 :       arg_ders[ ik ]=0;
     184           0 :       double dp_i=vals[ik]->difference( reference_args[i], arg[ik] );
     185           0 :       for(unsigned j=0; j<reference_args.size(); ++j) {
     186             :         double dp_j;
     187           0 :         unsigned jk=arg_der_index[j];
     188           0 :         if(i==j) {
     189             :           dp_j=dp_i;
     190             :         } else {
     191           0 :           dp_j=vals[jk]->difference( reference_args[j], arg[jk] );
     192             :         }
     193             : 
     194           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
     195           0 :         r+=dp_i*dp_j*metric(i,j);
     196             :       }
     197             :     }
     198             :   } else {
     199     1128420 :     for(unsigned i=0; i<reference_args.size(); ++i) {
     200      828660 :       unsigned ik=arg_der_index[i];
     201      828660 :       double dp_i=vals[ik]->difference( reference_args[i], arg[ik] );
     202      828660 :       r+=weights[i]*dp_i*dp_i;
     203      828660 :       arg_ders[ik]=2.0*weights[i]*dp_i;
     204             :     }
     205             :   }
     206      299760 :   if(!squared) {
     207         572 :     r=std::sqrt(r);
     208         572 :     double ir=1.0/(2.0*r);
     209        1716 :     for(unsigned i=0; i<arg_ders.size(); ++i) {
     210        1144 :       myder.setArgumentDerivatives( i, arg_ders[i]*ir );
     211             :     }
     212             :   } else {
     213     1126704 :     for(unsigned i=0; i<arg_ders.size(); ++i) {
     214             :       myder.setArgumentDerivatives( i, arg_ders[i] );
     215             :     }
     216             :   }
     217      299760 :   return r;
     218             : }
     219             : 
     220        1964 : void ReferenceArguments::extractArgumentDisplacement( const std::vector<Value*>& vals, const std::vector<double>& arg, std::vector<double>& dirout ) const {
     221        1964 :   if( hasmetric ) {
     222           0 :     plumed_error();
     223             :   } else {
     224        5892 :     for(unsigned j=0; j<reference_args.size(); ++j) {
     225        3928 :       unsigned jk=arg_der_index[j];
     226        3928 :       dirout[jk]=sqrtweight[j]*vals[jk]->difference( reference_args[j], arg[jk] );
     227             :     }
     228             :   }
     229        1964 : }
     230             : 
     231        2386 : double ReferenceArguments::projectArgDisplacementOnVector( const std::vector<double>& eigv, const std::vector<Value*>& vals, const std::vector<double>& arg, ReferenceValuePack& mypack ) const {
     232        2386 :   if( hasmetric ) {
     233           0 :     plumed_error();
     234             :   } else {
     235             :     double proj=0;
     236        7158 :     for(unsigned j=0; j<reference_args.size(); ++j) {
     237        4772 :       unsigned jk=arg_der_index[j];
     238        4772 :       proj += eigv[j]*sqrtweight[j]*vals[jk]->difference( reference_args[j], arg[jk] );
     239        4772 :       mypack.setArgumentDerivatives( jk, eigv[j]*sqrtweight[j] );
     240             :     }
     241        2386 :     return proj;
     242             :   }
     243             : }
     244             : 
     245         500 : void ReferenceArguments::displaceReferenceArguments( const double& weight, const std::vector<double>& displace ) {
     246             :   plumed_dbg_assert( displace.size()==getNumberOfReferenceArguments() );
     247        1466 :   for(unsigned i=0; i<displace.size(); ++i) {
     248         966 :     reference_args[i] += weight*displace[i];
     249             :   }
     250         500 : }
     251             : 
     252             : }