Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2019 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 "core/ActionRegister.h"
23 : #include "tools/Pbc.h"
24 : #include "tools/SwitchingFunction.h"
25 : #include "ActionVolume.h"
26 :
27 : //+PLUMEDOC VOLUMES INSPHERE
28 : /*
29 : This quantity can be used to calculate functions of the distribution of collective
30 : variables for the atoms that lie in a particular, user-specified part of of the cell.
31 :
32 : Each of the base quantities calculated by a multicolvar can can be assigned to a particular point in three
33 : dimensional space. For example, if we have the coordination numbers for all the atoms in the
34 : system each coordination number can be assumed to lie on the position of the central atom.
35 : Because each base quantity can be assigned to a particular point in space we can calculate functions of the
36 : distribution of base quantities in a particular part of the box by using:
37 :
38 : \f[
39 : \overline{s}_{\tau} = \frac{ \sum_i f(s_i) \sigma(r) }{ \sum_i \sigma(r) }
40 : \f]
41 :
42 : where the sum is over the collective variables, \f$s_i\f$, each of which can be thought to be at \f$ (x_i,y_i,z_i)\f$.
43 : The function \f$\sigma\f$ is a \ref switchingfunction that acts on the distance between the point at which the
44 : collective is located \f$(x_i,y_i,z_i)\f$ and the position of the atom that was specified using the ORIGIN keyword.
45 : In other words:
46 : \f[
47 : r = sqrt{ ( x_i - x_0)^2 + ( y_i - y_0)^2 + ( z_i - z_0)^2}
48 : \f]
49 : In short this function, \f$\sigma(r_{xy})\f$, measures whether or not the CV is within a sphere that is
50 : centered on the position of the atom specified using the keyword ORIGIN.
51 :
52 : The function \f$(s_i)\f$ can be any of the usual LESS_THAN, MORE_THAN, WITHIN etc that are used in all other multicolvars.
53 :
54 : When INCYLINDER is used with the \ref DENSITY action the number of atoms in the specified region is calculated
55 :
56 : \par Examples
57 :
58 : The input below can be use to calculate the average coordination numbers for those atoms that are within a sphere
59 : of radius 1.5 nm that is centered on the position of atom 101.
60 :
61 : \plumedfile
62 : c1: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=0.1}
63 : d2: INSPHERE ATOM=101 DATA=d1 RADIUS={TANH R_0=1.5} SIGMA=0.1 LOWER=-0.1 UPPER=0.1 MEAN
64 : PRINT ARG=d2.* FILE=colvar
65 : \endplumedfile
66 :
67 : */
68 : //+ENDPLUMEDOC
69 :
70 : namespace PLMD {
71 : namespace multicolvar {
72 :
73 6 : class VolumeInSphere : public ActionVolume {
74 : private:
75 : Vector origin;
76 : SwitchingFunction switchingFunction;
77 : public:
78 : static void registerKeywords( Keywords& keys );
79 : explicit VolumeInSphere(const ActionOptions& ao);
80 : void setupRegions();
81 : double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const ;
82 : };
83 :
84 6455 : PLUMED_REGISTER_ACTION(VolumeInSphere,"INSPHERE")
85 :
86 4 : void VolumeInSphere::registerKeywords( Keywords& keys ) {
87 4 : ActionVolume::registerKeywords( keys );
88 16 : keys.add("atoms","ATOM","the atom whose vicinity we are interested in examining");
89 16 : keys.add("compulsory","RADIUS","the switching function that tells us the extent of the sphereical region of interest");
90 8 : keys.remove("SIGMA");
91 4 : }
92 :
93 3 : VolumeInSphere::VolumeInSphere(const ActionOptions& ao):
94 : Action(ao),
95 3 : ActionVolume(ao)
96 : {
97 : std::vector<AtomNumber> atom;
98 6 : parseAtomList("ATOM",atom);
99 3 : if( atom.size()!=1 ) error("should only be one atom specified");
100 6 : log.printf(" center of sphere is at position of atom : %d\n",atom[0].serial() );
101 :
102 6 : std::string sw, errors; parse("RADIUS",sw);
103 3 : if(sw.length()==0) error("missing RADIUS keyword");
104 3 : switchingFunction.set(sw,errors);
105 3 : if( errors.length()!=0 ) error("problem reading RADIUS keyword : " + errors );
106 9 : log.printf(" radius of sphere is given by %s \n", ( switchingFunction.description() ).c_str() );
107 :
108 3 : checkRead(); requestAtoms(atom);
109 3 : }
110 :
111 82 : void VolumeInSphere::setupRegions() { }
112 :
113 155474 : double VolumeInSphere::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const {
114 : // Calculate position of atom wrt to origin
115 310948 : Vector fpos=pbcDistance( getPosition(0), cpos );
116 155474 : double dfunc, value = switchingFunction.calculateSqr( fpos.modulo2(), dfunc );
117 310948 : derivatives.zero(); derivatives = dfunc*fpos; refders[0] = -derivatives;
118 : // Add a virial contribution
119 155474 : vir -= Tensor(fpos,derivatives);
120 155474 : return value;
121 : }
122 :
123 : }
124 4839 : }
|
