Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2018 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 "MultiColvar.h"
23 : #include "tools/SwitchingFunction.h"
24 : #include "core/ActionRegister.h"
25 : #include "vesselbase/LessThan.h"
26 : #include "vesselbase/Between.h"
27 : #include "tools/Angle.h"
28 :
29 : #include <string>
30 : #include <cmath>
31 :
32 : using namespace std;
33 :
34 : namespace PLMD {
35 : namespace multicolvar {
36 :
37 : //+PLUMEDOC MCOLVAR INPLANEDISTANCES
38 : /*
39 : Calculate distances in the plane perpendicular to an axis
40 :
41 : Each quantity calculated in this CV uses the positions of two atoms, this indices of which are specified using the VECTORSTART and VECTOREND keywords, to specify the
42 : orientation of a vector, \f$\mathbf{n}\f$. The perpendicular distance between this vector and the position of some third atom is then computed using:
43 : \f[
44 : x_j = |\mathbf{r}_{j}| \sin (\theta_j)
45 : \f]
46 : where \f$\mathbf{r}_j\f$ is the distance between one of the two atoms that define the vector \f$\mathbf{n}\f$ and a third atom (atom \f$j\f$) and where \f$\theta_j\f$
47 : is the angle between the vector \f$\mathbf{n}\f$ and the vector \f$\mathbf{r}_{j}\f$. The \f$x_j\f$ values for each of the atoms specified using the GROUP keyword are calculated.
48 : Keywords such as MORE_THAN and LESS_THAN can then be used to calculate the number of these quantities that are more or less than a given cutoff.
49 :
50 : \par Examples
51 :
52 : The following input can be used to calculate the number of atoms that have indices greater than 3 and less than 101 that
53 : are within a cylinder with a radius of 0.3 nm that has its long axis aligned with the vector connecting atoms 1 and 2.
54 :
55 : \verbatim
56 : d1: INPLANEDISTANCES VECTORSTART=1 VECTOREND=2 GROUP=3-100 LESS_THAN={RATIONAL D_0=0.2 R_0=0.1}
57 : PRINT ARG=d1.lessthan FILE=colvar
58 : \endverbatim
59 :
60 :
61 : */
62 : //+ENDPLUMEDOC
63 :
64 0 : class InPlaneDistances : public MultiColvar {
65 : public:
66 : static void registerKeywords( Keywords& keys );
67 : explicit InPlaneDistances(const ActionOptions&);
68 : // active methods:
69 : virtual double compute(const unsigned& tindex, AtomValuePack& myatoms ) const ;
70 0 : bool isPeriodic() { return false; }
71 : };
72 :
73 2523 : PLUMED_REGISTER_ACTION(InPlaneDistances,"INPLANEDISTANCES")
74 :
75 1 : void InPlaneDistances::registerKeywords( Keywords& keys ) {
76 1 : MultiColvar::registerKeywords( keys );
77 1 : keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST");
78 1 : keys.use("MEAN"); keys.use("MIN"); keys.use("MAX"); keys.use("LESS_THAN");
79 1 : keys.use("MORE_THAN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
80 1 : keys.add("atoms","VECTORSTART","The first atom position that is used to define the normal to the plane of interest");
81 1 : keys.add("atoms","VECTOREND","The second atom position that is used to defin the normal to the plane of interest");
82 1 : keys.add("atoms-2","GROUP","The set of atoms for which you wish to calculate the in plane distance ");
83 1 : }
84 :
85 0 : InPlaneDistances::InPlaneDistances(const ActionOptions&ao):
86 0 : PLUMED_MULTICOLVAR_INIT(ao)
87 : {
88 : // Read in the atoms
89 0 : std::vector<AtomNumber> all_atoms;
90 0 : readThreeGroups("GROUP","VECTORSTART","VECTOREND",false,false,all_atoms);
91 :
92 : // Check atoms are OK
93 0 : if( getFullNumberOfTasks()!=getNumberOfAtoms()-2 ) error("you should specify one atom for VECTORSTART and one atom for VECTOREND only");
94 :
95 : // Setup the multicolvar base
96 0 : setupMultiColvarBase( all_atoms ); readVesselKeywords();
97 : // And check everything has been read in correctly
98 0 : checkRead();
99 :
100 : // Now check if we can use link cells
101 0 : bool use_link=false; double rcut;
102 0 : if( getNumberOfVessels()>0 ) {
103 0 : vesselbase::LessThan* lt=dynamic_cast<vesselbase::LessThan*>( getPntrToVessel(0) );
104 0 : if( lt ) {
105 0 : use_link=true; rcut=lt->getCutoff();
106 : } else {
107 0 : vesselbase::Between* bt=dynamic_cast<vesselbase::Between*>( getPntrToVessel(0) );
108 0 : if( bt ) use_link=true; rcut=bt->getCutoff();
109 : }
110 0 : if( use_link ) {
111 0 : for(unsigned i=1; i<getNumberOfVessels(); ++i) {
112 0 : vesselbase::LessThan* lt2=dynamic_cast<vesselbase::LessThan*>( getPntrToVessel(i) );
113 0 : vesselbase::Between* bt=dynamic_cast<vesselbase::Between*>( getPntrToVessel(i) );
114 0 : if( lt2 ) {
115 0 : double tcut=lt2->getCutoff();
116 0 : if( tcut>rcut ) rcut=tcut;
117 0 : } else if( bt ) {
118 0 : double tcut=bt->getCutoff();
119 0 : if( tcut>rcut ) rcut=tcut;
120 : } else {
121 0 : use_link=false;
122 : }
123 : }
124 : }
125 0 : if( use_link ) setLinkCellCutoff( rcut );
126 0 : }
127 0 : }
128 :
129 0 : double InPlaneDistances::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
130 0 : Vector normal=getSeparation( myatoms.getPosition(1), myatoms.getPosition(2) );
131 0 : Vector dir=getSeparation( myatoms.getPosition(1), myatoms.getPosition(0) );
132 0 : PLMD::Angle a; Vector ddij, ddik; double angle=a.compute(normal,dir,ddij,ddik);
133 0 : double sangle=sin(angle), cangle=cos(angle);
134 0 : double dd=dir.modulo(), invdd=1.0/dd, val=dd*sangle;
135 :
136 0 : addAtomDerivatives( 1, 0, dd*cangle*ddik + sangle*invdd*dir, myatoms );
137 0 : addAtomDerivatives( 1, 1, -dd*cangle*(ddik+ddij) - sangle*invdd*dir, myatoms );
138 0 : addAtomDerivatives( 1, 2, dd*cangle*ddij, myatoms );
139 0 : myatoms.addBoxDerivatives( 1, -dd*cangle*(Tensor(normal,ddij)+Tensor(dir,ddik)) - sangle*invdd*Tensor(dir,dir) );
140 :
141 0 : return val;
142 : }
143 :
144 : }
145 2523 : }
|
