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