Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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 "HBPammObject.h"
      23             : #include "tools/IFile.h"
      24             : 
      25             : namespace PLMD {
      26             : namespace pamm {
      27             : 
      28           4 : void HBPammObject::setup( const std::string& filename, const double& reg, multicolvar::MultiColvarBase* mybase, std::string& errorstr ) {
      29           4 :   mymulti=mybase; std::vector<std::string> valnames(3);
      30             :   valnames[0]="ptc"; valnames[1]="ssc"; valnames[2]="adc";
      31           4 :   std::vector<std::string> min(3), max(3); std::vector<bool> pbc(3, false);
      32           4 :   mypamm.setup( filename, reg, valnames, pbc, min, max, errorstr );
      33           4 : }
      34             : 
      35           4 : double HBPammObject::get_cutoff() const {
      36             :   double sfmax=0;
      37          48 :   for(unsigned k=0; k<mypamm.getNumberOfKernels(); ++k) {
      38          88 :     double rcut = mypamm.getKernelCenter(k)[2] + mypamm.getKernelSupport(k)[2];
      39          44 :     if( rcut>sfmax ) { sfmax=rcut; }
      40             :   }
      41           4 :   return sfmax;
      42             : }
      43             : 
      44      787016 : double HBPammObject::evaluate( const unsigned& dno, const unsigned& ano, const unsigned& hno,
      45             :                                const Vector& d_da, const double& md_da, multicolvar::AtomValuePack& myatoms ) const {
      46      787016 :   Vector d_dh = mymulti->getSeparation( myatoms.getPosition(dno), myatoms.getPosition(hno) ); double md_dh = d_dh.modulo(); // hydrogen - donor
      47      787016 :   Vector d_ah = mymulti->getSeparation( myatoms.getPosition(ano), myatoms.getPosition(hno) ); double md_ah = d_ah.modulo(); // hydrogen - acceptor
      48             : 
      49             :   // Create some vectors locally for pamm evaluation
      50      787016 :   std::vector<double> invals( 3 ), outvals( mypamm.getNumberOfKernels() );
      51      787016 :   std::vector<std:: vector<double> > der( mypamm.getNumberOfKernels() );
      52     9444192 :   for(unsigned i=0; i<der.size(); ++i) der[i].resize(3);
      53             : 
      54             :   // Evaluate the pamm object
      55      787016 :   invals[0]=md_dh - md_ah; invals[1]=md_dh+md_ah; invals[2]=md_da;
      56      787016 :   mypamm.evaluate( invals, outvals, der );
      57             : 
      58      787016 :   if( !mymulti->doNotCalculateDerivatives() ) {
      59          72 :     mymulti->addAtomDerivatives( 1, dno, ((-der[0][0])/md_dh)*d_dh, myatoms );
      60          72 :     mymulti->addAtomDerivatives( 1, ano, ((+der[0][0])/md_ah)*d_ah, myatoms  );
      61          72 :     mymulti->addAtomDerivatives( 1, hno, ((+der[0][0])/md_dh)*d_dh - ((+der[0][0])/md_ah)*d_ah, myatoms );
      62          72 :     myatoms.addBoxDerivatives( 1, ((-der[0][0])/md_dh)*Tensor(d_dh,d_dh) - ((-der[0][0])/md_ah)*Tensor(d_ah,d_ah) );
      63          72 :     mymulti->addAtomDerivatives( 1, dno, ((-der[0][1])/md_dh)*d_dh, myatoms );
      64          72 :     mymulti->addAtomDerivatives( 1, ano, ((-der[0][1])/md_ah)*d_ah, myatoms );
      65          72 :     mymulti->addAtomDerivatives( 1, hno, ((+der[0][1])/md_dh)*d_dh + ((+der[0][1])/md_ah)*d_ah, myatoms );
      66          72 :     myatoms.addBoxDerivatives( 1, ((-der[0][1])/md_dh)*Tensor(d_dh,d_dh) + ((-der[0][1])/md_ah)*Tensor(d_ah,d_ah) );
      67          72 :     mymulti->addAtomDerivatives( 1, dno, ((-der[0][2])/md_da)*d_da, myatoms );
      68          72 :     mymulti->addAtomDerivatives( 1, ano, ((+der[0][2])/md_da)*d_da, myatoms );
      69          72 :     myatoms.addBoxDerivatives( 1, ((-der[0][2])/md_da)*Tensor(d_da,d_da) );
      70             :   }
      71      787016 :   return outvals[0];
      72             : 
      73      787016 : }
      74             : 
      75             : }
      76             : }