Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2020-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 "Colvar.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/Angle.h"
25 :
26 : namespace PLMD {
27 : namespace colvar {
28 :
29 : //+PLUMEDOC COLVAR PROJECTION_ON_AXIS
30 : /*
31 : Calculate a position based on the projection along and extension from a defined axis.
32 :
33 : This variable takes 3 input atoms or pseudoatoms. It uses the two AXIS_ATOMS to define a linear vector.
34 : The position of the ATOM is then calculated relative to this vector, with two output components.
35 : The projection on the axis (proj) is the distance along the axis from the ATOM to the origin.
36 : The extension (ext) is the orthogonal distance between the ATOM and the axis.
37 :
38 : ## Examples
39 :
40 : This command tells plumed to define an axis, by calculating a vector that passes through atom 1 and atom 2.
41 : The position of atom 3 as a projection along this vector is calculated and printed to COLVAR1.
42 : At the same time, the perpendicular distance of atom 3 from the axis, the extension, is printed to COLVAR2.
43 :
44 : ```plumed
45 : poa: PROJECTION_ON_AXIS AXIS_ATOMS=1,2 ATOM=3
46 : PRINT ARG=poa.proj FILE=COLVAR1
47 : PRINT ARG=poa.ext FILE=COLVAR2
48 : ```
49 :
50 : A particular application of this variable could be to study the motion of a ligand relative to its binding pocket on a protein.
51 : In this set of commands, the anchor points a1 and a2 are defined using example atom numbers within the protein.
52 : As a2 is attempting to be as close as possible to the center of the binding pocket, a COM is used when there are no suitable protein atoms.
53 : Similarly, a COM is used to define the position of the ligand in lig1.
54 : The calculated projection of lig1 along the axis defined between a1 and a2 is printed to COLVAR1.
55 : The calculated perpendicular extension of lig1 from the axis defined between a1 and a2 is printed to COLVAR2.
56 :
57 : ```plumed
58 : a1: GROUP ATOMS=3754 # Anchor point 1
59 : a2: COM ATOMS=3019,4329,4744 # Anchor point 2
60 : lig1: COM ATOMS=5147-5190 # Ligand
61 : pp: PROJECTION_ON_AXIS AXIS_ATOMS=a1,a2 ATOM=lig1
62 : PRINT ARG=pp.proj FILE=COLVAR1
63 : PRINT ARG=pp.ext FILE=COLVAR2
64 : ```
65 :
66 : */
67 : //+ENDPLUMEDOC
68 :
69 : class ProjectionOnAxis : public Colvar {
70 : bool pbc;
71 :
72 : public:
73 : explicit ProjectionOnAxis(const ActionOptions&);
74 : // Active methods:
75 : virtual void calculate();
76 : static void registerKeywords( Keywords& keys );
77 : };
78 :
79 : PLUMED_REGISTER_ACTION(ProjectionOnAxis,"PROJECTION_ON_AXIS")
80 :
81 3 : void ProjectionOnAxis::registerKeywords( Keywords& keys ) {
82 3 : Colvar::registerKeywords(keys);
83 3 : keys.add("atoms","AXIS_ATOMS","The atoms that define the direction of the axis of interest");
84 3 : keys.add("atoms","ATOM","The atom whose position we want to project on the axis of interest");
85 6 : keys.addOutputComponent("proj","COMPONENTS","scalar","The value of the projection along the axis");
86 6 : keys.addOutputComponent("ext","COMPONENTS","scalar","The value of the extension from the axis");
87 6 : keys.setValueDescription("scalar","the value of the projection along the axis");
88 3 : }
89 :
90 1 : ProjectionOnAxis::ProjectionOnAxis(const ActionOptions&ao):
91 : PLUMED_COLVAR_INIT(ao),
92 1 : pbc(true) {
93 : std::vector<AtomNumber> axis_atoms;
94 2 : parseAtomList("AXIS_ATOMS",axis_atoms);
95 1 : if( axis_atoms.size()!=2 ) {
96 0 : error("There should only be two atoms specified to AXIS_ATOMS keyword");
97 : }
98 : std::vector<AtomNumber> atom;
99 2 : parseAtomList("ATOM",atom);
100 1 : if( atom.size()!=1 ) {
101 0 : error("There should only be one atom specified to ATOM keyword");
102 : }
103 1 : log.printf(" calculating projection of vector connecting atom %d and atom %d on vector connecting atom %d and atom %d \n",
104 : axis_atoms[0].serial(), atom[0].serial(), axis_atoms[0].serial(), axis_atoms[1].serial() );
105 1 : bool nopbc=!pbc;
106 1 : parseFlag("NOPBC",nopbc);
107 1 : pbc=!nopbc;
108 :
109 1 : if(pbc) {
110 1 : log.printf(" using periodic boundary conditions\n");
111 : } else {
112 0 : log.printf(" not using periodic boundary conditions\n");
113 : }
114 :
115 : // Add values to store data
116 2 : addComponentWithDerivatives("proj");
117 2 : componentIsNotPeriodic("proj");
118 2 : addComponentWithDerivatives("ext");
119 2 : componentIsNotPeriodic("ext");
120 : // Get all the atom positions
121 1 : axis_atoms.push_back( atom[0] );
122 1 : requestAtoms(axis_atoms);
123 1 : checkRead();
124 1 : }
125 :
126 : // Calculator
127 25 : void ProjectionOnAxis::calculate() {
128 :
129 25 : Vector rik, rjk;
130 25 : if( pbc ) {
131 25 : rik = pbcDistance( getPosition(2), getPosition(0) );
132 25 : rjk = pbcDistance( getPosition(2), getPosition(1) );
133 : } else {
134 0 : rik = delta( getPosition(2), getPosition(0) );
135 0 : rjk = delta( getPosition(2), getPosition(1) );
136 : }
137 25 : Vector rij = delta( rik, rjk );
138 25 : double dij = rij.modulo();
139 25 : Vector nij = (1.0/dij)*rij;
140 25 : Tensor dij_a1;
141 : // Derivative of director connecting atom1 - atom2 wrt the position of atom 1
142 25 : dij_a1(0,0) = ( -(nij[1]*nij[1]+nij[2]*nij[2])/dij ); // dx/dx
143 25 : dij_a1(0,1) = ( nij[0]*nij[1]/dij ); // dx/dy
144 25 : dij_a1(0,2) = ( nij[0]*nij[2]/dij ); // dx/dz
145 25 : dij_a1(1,0) = ( nij[1]*nij[0]/dij ); // dy/dx
146 25 : dij_a1(1,1) = ( -(nij[0]*nij[0]+nij[2]*nij[2])/dij ); // dy/dy
147 25 : dij_a1(1,2) = ( nij[1]*nij[2]/dij );
148 25 : dij_a1(2,0) = ( nij[2]*nij[0]/dij );
149 25 : dij_a1(2,1) = ( nij[2]*nij[1]/dij );
150 25 : dij_a1(2,2) = ( -(nij[1]*nij[1]+nij[0]*nij[0])/dij );
151 :
152 : // Calculate dot product and derivatives
153 25 : double d = dotProduct( -rik, nij );
154 25 : Vector dd1 = matmul(-rik, dij_a1) - nij;
155 25 : Vector dd2 = matmul(rik, dij_a1);
156 25 : Vector dd3 = nij;
157 50 : Value* pval=getPntrToComponent("proj");
158 : pval->set( d );
159 25 : setAtomsDerivatives( pval, 0, dd1 );
160 25 : setAtomsDerivatives( pval, 1, dd2 );
161 25 : setAtomsDerivatives( pval, 2, dd3 );
162 25 : setBoxDerivatives( pval, -Tensor( rik, dd1 ) - Tensor( rjk, dd2 ) );
163 : // Calculate derivatives of perpendicular distance from axis
164 25 : double c = std::sqrt( rik.modulo2() - d*d );
165 25 : double invc = (1.0/c);
166 : // Calculate derivatives of the other thing
167 25 : Vector der1 = invc*(rik - d*dd1);
168 25 : Vector der2 = invc*(-d*dd2);
169 25 : Vector der3 = invc*(-rik - d*dd3);
170 :
171 50 : Value* cval=getPntrToComponent("ext");
172 : cval->set( c );
173 25 : setAtomsDerivatives( cval, 0, der1 );
174 25 : setAtomsDerivatives( cval, 1, der2 );
175 25 : setAtomsDerivatives( cval, 2, der3 );
176 25 : setBoxDerivatives( cval, -Tensor( rik, der1 ) - Tensor( rjk, der2 ) );
177 25 : }
178 :
179 : }
180 : }
181 :
182 :
183 :
|
