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