Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-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 "DistanceFromContourBase.h"
      23             : 
      24             : namespace PLMD {
      25             : namespace contour {
      26             : 
      27           6 : void DistanceFromContourBase::registerKeywords( Keywords& keys ) {
      28           6 :   Action::registerKeywords( keys ); ActionWithValue::registerKeywords( keys );
      29           6 :   ActionAtomistic::registerKeywords( keys ); ActionWithArguments::registerKeywords( keys );
      30          12 :   keys.remove("NUMERICAL_DERIVATIVES"); keys.use("ARG");
      31          12 :   keys.add("atoms","POSITIONS","the positions of the atoms that we are calculating the contour from");
      32          12 :   keys.add("atoms","ATOM","The atom whose perpendicular distance we are calculating from the contour");
      33          12 :   keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
      34          12 :   keys.add("compulsory","KERNEL","GAUSSIAN","the kernel function you are using.  More details on  the kernels available "
      35             :            "in plumed plumed can be found in \\ref kernelfunctions.");
      36          12 :   keys.add("compulsory","CUTOFF","6.25","the cutoff at which to stop evaluating the kernel functions is set equal to sqrt(2*x)*bandwidth in each direction where x is this number");
      37          12 :   keys.add("compulsory","CONTOUR","the value we would like for the contour");
      38           6 : }
      39             : 
      40           2 : DistanceFromContourBase::DistanceFromContourBase( const ActionOptions& ao ):
      41             :   Action(ao),
      42             :   ActionWithValue(ao),
      43             :   ActionAtomistic(ao),
      44             :   ActionWithArguments(ao),
      45             :   mymin(this),
      46           2 :   nactive(0)
      47             : {
      48           2 :   if( getNumberOfArguments()>1 ) error("should only use one argument for this action");
      49           2 :   if( getNumberOfArguments()==1 ) {
      50           1 :     if( getPntrToArgument(0)->getRank()!=1 ) error("ARG for distance from contour should be rank one");
      51           1 :     getPntrToArgument(0)->buildDataStore();
      52             :   }
      53             :   // Read in the multicolvar/atoms
      54           4 :   std::vector<AtomNumber> atoms; parseAtomList("POSITIONS",atoms);
      55           4 :   std::vector<AtomNumber> origin; parseAtomList("ATOM",origin);
      56           2 :   if( origin.size()!=1 ) error("should only specify one atom for origin keyword");
      57             :   std::vector<AtomNumber> center;
      58           4 :   if( keywords.exists("ORIGIN") ) {
      59           2 :     parseAtomList("ORIGIN",center);
      60           1 :     if( center.size()!=1 ) error("should only specify one atom for center keyword");
      61             :   }
      62             : 
      63           2 :   if( center.size()==1 ) log.printf("  calculating distance between atom %d and contour along vector connecting it to atom %d \n", origin[0].serial(),center[0].serial() );
      64           1 :   else log.printf("  calculating distance between atom %d and contour \n", origin[0].serial() );
      65             : 
      66           2 :   log.printf("  contour is in field constructed from positions of atoms : ");
      67         515 :   for(unsigned i=0; i<atoms.size(); ++i) log.printf("%d ",atoms[i].serial() );
      68           2 :   if( getNumberOfArguments()==1 ) {
      69           1 :     if( getPntrToArgument(0)->getShape()[0]!=atoms.size() ) error("mismatch between number of atoms and size of vector specified using ARG keyword");
      70           1 :     log.printf("\n  and weights from %s \n", getPntrToArgument(0)->getName().c_str() );
      71             :   } else {
      72           1 :     log.printf("\n  all weights are set equal to one \n");
      73             :   }
      74             :   // Request everything we need
      75           2 :   active_list.resize( atoms.size(), 0 );
      76           2 :   std::vector<Value*> args( getArguments() ); atoms.push_back( origin[0] );
      77           2 :   if( center.size()==1 ) atoms.push_back( center[0] );
      78           2 :   requestArguments( args ); requestAtoms( atoms );
      79             :   // Fix to request arguments
      80           2 :   if( args.size()==1 ) addDependency( args[0]->getPntrToAction() );
      81             : 
      82             :   // Read in details of phase field construction
      83           8 :   parseVector("BANDWIDTH",bw); parse("KERNEL",kerneltype); parse("CONTOUR",contour);
      84           4 :   std::string errors; switchingFunction.set( kerneltype + " R_0=1.0 NOSTRETCH", errors );
      85           2 :   if( errors.length()!=0 ) error("problem reading switching function description " + errors);
      86           8 :   double det=1; for(unsigned i=0; i<bw.size(); ++i) det*=bw[i]*bw[i];
      87           2 :   gvol=1.0; if( kerneltype=="GAUSSIAN" ) gvol=pow( 2*pi, 0.5*bw.size() ) * pow( det, 0.5 );
      88           2 :   log.printf("  constructing phase field using %s kernels with bandwidth (%f, %f, %f) \n",kerneltype.c_str(), bw[0], bw[1], bw[2] );
      89             : 
      90             :   // And a cutoff
      91           2 :   std::vector<double> pp( bw.size(),0 );
      92           4 :   double dp2cutoff; parse("CUTOFF",dp2cutoff); double rcut =  sqrt(2*dp2cutoff)*bw[0];
      93           6 :   for(unsigned j=1; j<bw.size(); ++j) {
      94           4 :     if( sqrt(2*dp2cutoff)*bw[j]>rcut ) rcut=sqrt(2*dp2cutoff)*bw[j];
      95             :   }
      96           2 :   rcut2=rcut*rcut;
      97             : 
      98             :   // Create the vector of values that holds the position
      99           2 :   forcesToApply.resize( 3*getNumberOfAtoms() + 9 );
     100           2 : }
     101             : 
     102         154 : void DistanceFromContourBase::lockRequests() {
     103             :   ActionWithArguments::lockRequests();
     104             :   ActionAtomistic::lockRequests();
     105         154 : }
     106             : 
     107         154 : void DistanceFromContourBase::unlockRequests() {
     108             :   ActionWithArguments::unlockRequests();
     109             :   ActionAtomistic::unlockRequests();
     110         154 : }
     111             : 
     112      190081 : double DistanceFromContourBase::evaluateKernel( const Vector& cpos, const Vector& apos, std::vector<double>& der ) const {
     113      190081 :   Vector distance = pbcDistance( getPosition(getNumberOfAtoms()-1), cpos );
     114             :   Vector dist2 = pbcDistance( distance, apos );
     115      760324 :   double dval=0; for(unsigned j=0; j<3; ++j) { der[j] = dist2[j]/bw[j]; dval += der[j]*der[j]; }
     116      190081 :   double dfunc, newval = switchingFunction.calculateSqr( dval, dfunc ) / gvol;
     117      760324 :   double tmp = dfunc / gvol; for(unsigned j=0; j<3; ++j) der[j] *= -tmp;
     118      190081 :   return newval;
     119             : }
     120             : 
     121        3283 : double DistanceFromContourBase::getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der ) {
     122             :   // Transer the position to the local Vector
     123       13132 :   for(unsigned j=0; j<3; ++j) pval[j] = x[j];
     124             :   // Now find the contour
     125        3283 :   double sumk = 0, sumd = 0; std::vector<double> pp(3), ddd(3,0);
     126      193090 :   for(unsigned i=0; i<nactive; ++i) {
     127      189807 :     double newval = evaluateKernel( getPosition(active_list[i]), pval, ddd );
     128             : 
     129      189807 :     if( getNumberOfArguments()==1 ) {
     130      186966 :       sumk += getPntrToArgument(0)->get(active_list[i])*newval;
     131      186966 :       sumd += newval;
     132        2841 :     } else sumk += newval;
     133             :   }
     134        3283 :   if( getNumberOfArguments()==0 ) return sumk - contour;
     135         442 :   if( fabs(sumk)<epsilon ) return -contour;
     136         195 :   return (sumk/sumd) - contour;
     137             : }
     138             : 
     139         154 : void DistanceFromContourBase::apply() {
     140         154 :   if( !checkForForces() ) return ;
     141         137 :   unsigned ind=0; setForcesOnAtoms( getForcesToApply(), ind );
     142             : }
     143             : 
     144             : }
     145             : }