Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 : Copyright (c) 2015-2020 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 : #include "VolumeShortcut.h" 27 : 28 : //+PLUMEDOC VOLUMES INSPHERE 29 : /* 30 : This quantity can be used to calculate functions of the distribution of collective 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=c1 RADIUS={TANH R_0=1.5} MEAN 64 : PRINT ARG=d2.* FILE=colvar 65 : \endplumedfile 66 : 67 : */ 68 : //+ENDPLUMEDOC 69 : 70 : //+PLUMEDOC MCOLVAR INSPHERE_CALC 71 : /* 72 : Calculate a vector from the input positions with elements equal to one when the positions are in a particular part of the cell and elements equal to zero otherwise 73 : 74 : \par Examples 75 : 76 : */ 77 : //+ENDPLUMEDOC 78 : 79 : namespace PLMD { 80 : namespace volumes { 81 : 82 : class VolumeInSphere : public ActionVolume { 83 : private: 84 : Vector origin; 85 : SwitchingFunction switchingFunction; 86 : public: 87 : static void registerKeywords( Keywords& keys ); 88 : explicit VolumeInSphere(const ActionOptions& ao); 89 : void setupRegions() override; 90 : double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override; 91 : }; 92 : 93 : PLUMED_REGISTER_ACTION(VolumeInSphere,"INSPHERE_CALC") 94 : char glob_sphere[] = "INSPHERE"; 95 : typedef VolumeShortcut<glob_sphere> VolumeInSphereShortcut; 96 : PLUMED_REGISTER_ACTION(VolumeInSphereShortcut,"INSPHERE") 97 : 98 36 : void VolumeInSphere::registerKeywords( Keywords& keys ) { 99 36 : ActionVolume::registerKeywords( keys ); keys.setDisplayName("INSPHERE"); 100 72 : keys.add("atoms","CENTER","the atom whose vicinity we are interested in examining"); 101 72 : keys.add("atoms-2","ATOM","the atom whose vicinity we are interested in examining"); 102 72 : keys.add("compulsory","RADIUS","the switching function that tells us the extent of the sphereical region of interest"); 103 36 : keys.remove("SIGMA"); 104 36 : } 105 : 106 10 : VolumeInSphere::VolumeInSphere(const ActionOptions& ao): 107 : Action(ao), 108 10 : ActionVolume(ao) 109 : { 110 20 : std::vector<AtomNumber> atom; parseAtomList("CENTER",atom); 111 10 : if( atom.size()==0 ) parseAtomList("ATOM",atom); 112 10 : if( atom.size()!=1 ) error("should only be one atom specified"); 113 10 : log.printf(" center of sphere is at position of atom : %d\n",atom[0].serial() ); 114 : 115 20 : std::string sw, errors; parse("RADIUS",sw); 116 10 : if(sw.length()==0) error("missing RADIUS keyword"); 117 10 : switchingFunction.set(sw,errors); 118 10 : if( errors.length()!=0 ) error("problem reading RADIUS keyword : " + errors ); 119 10 : log.printf(" radius of sphere is given by %s \n", ( switchingFunction.description() ).c_str() ); 120 : 121 10 : checkRead(); requestAtoms(atom); 122 10 : } 123 : 124 120 : void VolumeInSphere::setupRegions() { } 125 : 126 423194 : double VolumeInSphere::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const { 127 : // Calculate position of atom wrt to origin 128 423194 : Vector fpos=pbcDistance( getPosition(0), cpos ); 129 423194 : double dfunc, value = switchingFunction.calculateSqr( fpos.modulo2(), dfunc ); 130 423194 : derivatives.zero(); derivatives = dfunc*fpos; refders[0] = -derivatives; 131 : // Add a virial contribution 132 423194 : vir -= Tensor(fpos,derivatives); 133 423194 : return value; 134 : } 135 : 136 : } 137 : }