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