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             : #include "core/ActionRegister.h"
      24             : 
      25             : //+PLUMEDOC COLVAR DISTANCE_FROM_SPHERICAL_CONTOUR
      26             : /*
      27             : Calculate the perpendicular distance from a Willard-Chandler dividing surface.
      28             : 
      29             : This action works similarly to [DISTANCE_FROM_CONTOUR](DISTANCE_FROM_CONTOUR.md). Within this action a field is constructed that measures the density
      30             : of the system at each point in space using:
      31             : 
      32             : $$
      33             : p(x,y,x) = \sum_{i=1}^N K\left[\frac{x-x_i}{\sigma_x},\frac{y-y_i}{\sigma_y},\frac{z-z_i}{\sigma_z} \right]
      34             : $$
      35             : 
      36             : In this expression $\sigma_x, \sigma_y$ and $\sigma_z$ are bandwidth parameters and
      37             : $K$ is one of a Gaussian kernel function.  With that field in place we can define a Willard-Chandler
      38             : surface is defined a surface of constant density in the above field $p(x,y,z)$.
      39             : In other words, we can define a set of points, $(x',y',z')$, in the box which have:
      40             : 
      41             : $$
      42             : p(x',y',z') = \rho
      43             : $$
      44             : 
      45             : where $\rho$ is some target density. In [DISTANCE_FROM_CONTOUR](DISTANCE_FROM_CONTOUR.md) we assume that this set of points lie on a manifold that
      46             : has the same topology as one or multiple planes.  Here, by contrast, we assume that this set of points lie on a manifold that has the same topology
      47             : as a sphere.  This action then returns the distance between this spherical manifold and the position of a test particle.  This distance is measured
      48             : along a vector perpendicular to the manifold.
      49             : 
      50             : ## Examples
      51             : 
      52             : The following input calculates a [CONTACT_MATRIX](CONTACT_MATRIX.md) between a set of atoms in which element $i,j$ is only non-zero if atoms $i$ and $j$
      53             : are within 6 nm of each other.  We then use this matrix as input for a [DFSCLUSTERING](DFSCLUSTERING.md) action that finds the largest connected component
      54             : in the matrix.  The [CENTER](CENTER.md) of this cluster is then identified and the location of an isocontour in:
      55             : 
      56             : $$
      57             : p(x,y,x) = \sum_{i=1}^N \xi_i f(c_i) K\left[\frac{x-x_i}{\sigma_x},\frac{y-y_i}{\sigma_y},\frac{z-z_i}{\sigma_z} \right]
      58             : $$
      59             : 
      60             : is found using this action.  In this expression $\xi_i$ is 1 if atom $i$ is part of the largest cluster and zero otherwise, $c_i$ is the coordination number of atom $i$ and
      61             : $f$ is a swtiching function.  The distance between this isocontour and position of atom 513 as well as the distance between `com` (the center of the largest cluster)
      62             : and the isocontour is then output to a file called colvar.
      63             : 
      64             : ```plumed
      65             : ones: ONES SIZE=512
      66             : # Calculate contact matrix
      67             : c1_mat: CONTACT_MATRIX GROUP=1-512 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
      68             : # Calculate coordination numbers
      69             : c1: MATRIX_VECTOR_PRODUCT ARG=c1_mat,ones
      70             : # Select coordination numbers that are more than 2.0
      71             : cf: MORE_THAN ARG=c1 SWITCH={RATIONAL D_0=2.0 R_0=0.1}
      72             : # Find largest cluster
      73             : dfs: DFSCLUSTERING ARG=c1_mat
      74             : clust1: CLUSTER_WEIGHTS CLUSTERS=dfs CLUSTER=1
      75             : com: CENTER ATOMS=1-512 WEIGHTS=clust1 PHASES
      76             : # Filtered coordination numbers for atoms in largest cluster
      77             : ff: CUSTOM ARG=clust1,cf FUNC=x*y PERIODIC=NO
      78             : # Now do the multicolvar surface
      79             : dd: DISTANCE_FROM_SPHERICAL_CONTOUR ARG=ff POSITIONS=1-512 ATOM=513 ORIGIN=com BANDWIDTH=1.0,1.0,1.0 CONTOUR=0.5
      80             : PRINT ARG=dd.* FILE=colvar
      81             : ```
      82             : 
      83             : */
      84             : //+ENDPLUMEDOC
      85             : 
      86             : namespace PLMD {
      87             : namespace contour {
      88             : 
      89             : class DistanceFromSphericalContour : public DistanceFromContourBase {
      90             : public:
      91             :   static void registerKeywords( Keywords& keys );
      92             :   explicit DistanceFromSphericalContour( const ActionOptions& );
      93             :   void calculate();
      94             :   void evaluateDerivatives( const Vector& root1, const double& root2 );
      95             : };
      96             : 
      97             : PLUMED_REGISTER_ACTION(DistanceFromSphericalContour,"DISTANCE_FROM_SPHERICAL_CONTOUR")
      98             : 
      99           3 : void DistanceFromSphericalContour::registerKeywords( Keywords& keys ) {
     100           3 :   DistanceFromContourBase::registerKeywords( keys );
     101           6 :   keys.addOutputComponent("dist","default","scalar","the distance between the reference atom and the nearest contour");
     102           6 :   keys.addOutputComponent("radius","default","scalar","the radial distance from the center of the contour to the edge");
     103           3 :   keys.add("atoms","ORIGIN","The position of the center of the region that the contour encloses");
     104           3 :   keys.addDOI("10.1021/acs.jpcb.8b03661");
     105           3 : }
     106             : 
     107           1 : DistanceFromSphericalContour::DistanceFromSphericalContour( const ActionOptions& ao ):
     108             :   Action(ao),
     109           1 :   DistanceFromContourBase(ao) {
     110             :   // Create the values
     111             :   std::vector<std::size_t> shape;
     112           1 :   addComponentWithDerivatives("dist", shape );
     113           1 :   componentIsNotPeriodic("dist");
     114           1 :   addComponent("radius", shape );
     115           2 :   componentIsNotPeriodic("radius");
     116           1 : }
     117             : 
     118          17 : void DistanceFromSphericalContour::calculate() {
     119             :   // Check box is orthorhombic
     120          17 :   if( !getPbc().isOrthorombic() ) {
     121           0 :     error("cell box must be orthorhombic");
     122             :   }
     123             : 
     124             :   // Calculate the director of the vector connecting the center of the sphere to the molecule of interest
     125          17 :   Vector dirv = pbcDistance( getPosition(getNumberOfAtoms()-1), getPosition(getNumberOfAtoms()-2) );
     126          17 :   double len=dirv.modulo();
     127             :   dirv /= len;
     128             :   // Now work out which atoms need to be considered explicitly
     129          17 :   pbcDistance( getPosition(getNumberOfAtoms()-1), getPosition(0) );
     130          17 :   nactive=1;
     131          17 :   active_list[0]=0;
     132        8704 :   for(unsigned j=1; j<getNumberOfAtoms()-2; ++j) {
     133        8687 :     if( getNumberOfArguments()==1 ) {
     134        8687 :       if( getPntrToArgument(0)->get(j)<epsilon ) {
     135             :         continue;
     136             :       }
     137             :     }
     138        6698 :     active_list[nactive]=j;
     139        6698 :     nactive++;
     140        6698 :     Vector distance=pbcDistance( getPosition(getNumberOfAtoms()-1), getPosition(j) );
     141             :     double dp = dotProduct( distance, dirv );
     142        6698 :     double cp = distance.modulo2() - dp*dp;
     143        6698 :     if( cp<rcut2 ) {
     144         476 :       active_list[nactive]=j;
     145         476 :       nactive++;
     146             :     }
     147             :   }
     148             :   // Get maximum length to fit in box
     149          17 :   double hbox = 0.5*getBox()(0,0);
     150          17 :   if( 0.5*getBox()(1,1)<hbox ) {
     151           0 :     hbox = 0.5*getBox()(1,1);
     152             :   }
     153          17 :   if( 0.5*getBox()(2,2)<hbox ) {
     154           0 :     hbox = 0.5*getBox()(2,2);
     155             :   }
     156             :   // Set initial guess for position of contour to position of closest molecule in region
     157          17 :   std::vector<double> pos1(3), dirv2(3);
     158          68 :   for(unsigned k=0; k<3; ++k) {
     159          51 :     dirv2[k]=hbox*dirv[k];
     160          51 :     pos1[k]=0;
     161             :   }
     162             :   // Now do a search for the contours
     163             :   findContour( dirv2, pos1 );
     164             :   // Now find the distance between the center of the sphere and the contour
     165          17 :   double rad = sqrt( pos1[0]*pos1[0] + pos1[1]*pos1[1] + pos1[2]*pos1[2] );
     166             :   // Set the radius
     167          17 :   getPntrToComponent("radius")->set( rad );
     168             :   // Set the distance between the contour and the molecule
     169          17 :   getPntrToComponent("dist")->set( len - rad );
     170             : 
     171             :   // Now calculate the derivatives
     172          17 :   if( !doNotCalculateDerivatives() ) {
     173           0 :     plumed_merror("derivatives not implemented");
     174             :   }
     175          17 : }
     176             : 
     177           0 : void DistanceFromSphericalContour::evaluateDerivatives( const Vector& root1, const double& root2 ) {
     178           0 :   plumed_error();
     179             : }
     180             : 
     181             : }
     182             : }