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, using 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 : \par 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 : \plumedfile
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 : \endplumedfile
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 : \plumedfile
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 : \endplumedfile
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 6 : keys.add("atoms","AXIS_ATOMS","The atoms that define the direction of the axis of interest");
84 6 : 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 : {
94 : std::vector<AtomNumber> axis_atoms;
95 2 : parseAtomList("AXIS_ATOMS",axis_atoms);
96 1 : if( axis_atoms.size()!=2 ) error("There should only be two atoms specified to AXIS_ATOMS keyword");
97 : std::vector<AtomNumber> atom;
98 2 : parseAtomList("ATOM",atom);
99 1 : if( atom.size()!=1 ) error("There should only be one atom specified to ATOM keyword");
100 1 : log.printf(" calculating projection of vector connecting atom %d and atom %d on vector connecting atom %d and atom %d \n",
101 : axis_atoms[0].serial(), atom[0].serial(), axis_atoms[0].serial(), axis_atoms[1].serial() );
102 1 : bool nopbc=!pbc;
103 1 : parseFlag("NOPBC",nopbc);
104 1 : pbc=!nopbc;
105 :
106 1 : if(pbc) log.printf(" using periodic boundary conditions\n");
107 0 : else log.printf(" not using periodic boundary conditions\n");
108 :
109 : // Add values to store data
110 3 : addComponentWithDerivatives("proj"); componentIsNotPeriodic("proj");
111 3 : addComponentWithDerivatives("ext"); componentIsNotPeriodic("ext");
112 : // Get all the atom positions
113 1 : axis_atoms.push_back( atom[0] );
114 1 : requestAtoms(axis_atoms);
115 1 : checkRead();
116 1 : }
117 :
118 : // Calculator
119 25 : void ProjectionOnAxis::calculate() {
120 :
121 25 : Vector rik, rjk;
122 25 : if( pbc ) {
123 25 : rik = pbcDistance( getPosition(2), getPosition(0) );
124 25 : rjk = pbcDistance( getPosition(2), getPosition(1) );
125 : } else {
126 0 : rik = delta( getPosition(2), getPosition(0) );
127 0 : rjk = delta( getPosition(2), getPosition(1) );
128 : }
129 25 : Vector rij = delta( rik, rjk ); double dij = rij.modulo();
130 25 : Vector nij = (1.0/dij)*rij; Tensor dij_a1;
131 : // Derivative of director connecting atom1 - atom2 wrt the position of atom 1
132 25 : dij_a1(0,0) = ( -(nij[1]*nij[1]+nij[2]*nij[2])/dij ); // dx/dx
133 25 : dij_a1(0,1) = ( nij[0]*nij[1]/dij ); // dx/dy
134 25 : dij_a1(0,2) = ( nij[0]*nij[2]/dij ); // dx/dz
135 25 : dij_a1(1,0) = ( nij[1]*nij[0]/dij ); // dy/dx
136 25 : dij_a1(1,1) = ( -(nij[0]*nij[0]+nij[2]*nij[2])/dij ); // dy/dy
137 25 : dij_a1(1,2) = ( nij[1]*nij[2]/dij );
138 25 : dij_a1(2,0) = ( nij[2]*nij[0]/dij );
139 25 : dij_a1(2,1) = ( nij[2]*nij[1]/dij );
140 25 : dij_a1(2,2) = ( -(nij[1]*nij[1]+nij[0]*nij[0])/dij );
141 :
142 : // Calculate dot product and derivatives
143 25 : double d = dotProduct( -rik, nij );
144 25 : Vector dd1 = matmul(-rik, dij_a1) - nij;
145 25 : Vector dd2 = matmul(rik, dij_a1);
146 25 : Vector dd3 = nij;
147 50 : Value* pval=getPntrToComponent("proj"); pval->set( d );
148 25 : setAtomsDerivatives( pval, 0, dd1 );
149 25 : setAtomsDerivatives( pval, 1, dd2 );
150 25 : setAtomsDerivatives( pval, 2, dd3 );
151 25 : setBoxDerivatives( pval, -Tensor( rik, dd1 ) - Tensor( rjk, dd2 ) );
152 : // Calculate derivatives of perpendicular distance from axis
153 25 : double c = std::sqrt( rik.modulo2() - d*d ); double invc = (1.0/c);
154 : // Calculate derivatives of the other thing
155 25 : Vector der1 = invc*(rik - d*dd1);
156 25 : Vector der2 = invc*(-d*dd2);
157 25 : Vector der3 = invc*(-rik - d*dd3);
158 :
159 50 : Value* cval=getPntrToComponent("ext"); cval->set( c );
160 25 : setAtomsDerivatives( cval, 0, der1 );
161 25 : setAtomsDerivatives( cval, 1, der2 );
162 25 : setAtomsDerivatives( cval, 2, der3 );
163 25 : setBoxDerivatives( cval, -Tensor( rik, der1 ) - Tensor( rjk, der2 ) );
164 25 : }
165 :
166 : }
167 : }
168 :
169 :
170 :
|
