Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 "core/ActionWithMatrix.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/Torsion.h"
25 :
26 : //+PLUMEDOC MCOLVAR TORSIONS_MATRIX
27 : /*
28 : Calculate the matrix of torsions between two vectors of molecules
29 :
30 : This action was implemented to ensure that we can calculate the [SMAC](SMAC.md) collective variable that is discussed in
31 : [this paper](https://www.sciencedirect.com/science/article/abs/pii/S0009250914004503?via%3Dihub). This particular action
32 : tracks the relative orientations for all the pairs of molecules in a set much like the variables described in the crystdistrib module.
33 :
34 : The orientations of molecules can be specified using either [PLANE](PLANE.md) or [DISTANCE](DISTANCE.md). The example below shows how you can use
35 : internal vectors connecting two atoms in the molecules to define the orientation of that molecule. Three of these internal
36 : vectors are calculated using a DISTANCE command in the input below. The matrix of torsional angles between these various
37 : vectors is then computed:
38 :
39 : ```plumed
40 : d1: DISTANCE ATOMS1=1,5 ATOMS2=11,15 ATOMS3=21,25 COMPONENTS
41 : s: VSTACK ARG=d1.x,d1.y,d1.z
42 : sT: TRANSPOSE ARG=s
43 : m: TORSIONS_MATRIX ARG=s,sT POSITIONS1=1,11,21 POSITIONS2=1,11,21
44 : PRINT ARG=m FILE=matrix
45 : ```
46 :
47 : In this example, the torsional angle in element $(1,2)$ of the matrix with label `m` is the angle between the plane containing atoms 1,5 and 10 and the plane
48 : connecting atoms 1,10 and 15. In other words, the elements in this matrix are the torsional angles between the vectors in the input matrices
49 : around the vector connecting the corresponding atomic positions that are specified using the `POSTIONS` keyword.
50 :
51 : You can also calculate a matrix of torsional angles between two different groups of molecules by using an input like the one below:
52 :
53 : ```plumed
54 : pA: PLANE ATOMS1=1,2,3 ATOMS2=11,12,13
55 : sA: VSTACK ARG=pA.x,pA.y,pA.z
56 : pB: PLANE ATOMS1=21,22,23 ATOMS2=31,32,33 ATOMS3=41,42,43
57 : sB: VSTACK ARG=pB.x,pB.y,pB.z
58 : sBT: TRANSPOSE ARG=sB
59 : m: TORSIONS_MATRIX ARG=sA,sBT POSITIONS1=1,11 POSITIONS2=21,31,41
60 : PRINT ARG=m FILE=matrix
61 : ```
62 :
63 : In this example, the orientations of the molecules are specified using the [PLANE](PLANE.md) action and is given by a normal to the plane containing the three atoms from the molecule
64 : that was specified. The final output is $2 \times 3$ matrix that contains all the torsional angles between the molecules defined by the two PLANE actions.
65 :
66 : */
67 : //+ENDPLUMEDOC
68 :
69 : namespace PLMD {
70 : namespace adjmat {
71 :
72 : class TorsionsMatrix : public ActionWithMatrix {
73 : private:
74 : unsigned nderivatives;
75 : bool stored_matrix1, stored_matrix2;
76 : public:
77 : static void registerKeywords( Keywords& keys );
78 : explicit TorsionsMatrix(const ActionOptions&);
79 : unsigned getNumberOfDerivatives();
80 0 : unsigned getNumberOfColumns() const override {
81 0 : return getConstPntrToComponent(0)->getShape()[1];
82 : }
83 : void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
84 : void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
85 : void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
86 : };
87 :
88 : PLUMED_REGISTER_ACTION(TorsionsMatrix,"TORSIONS_MATRIX")
89 :
90 14 : void TorsionsMatrix::registerKeywords( Keywords& keys ) {
91 14 : ActionWithMatrix::registerKeywords(keys);
92 28 : keys.addInputKeyword("compulsory","ARG","matrix","an Nx3 and a 3xN matrix that contain the bond vectors that you would like to determine the torsion angles between");
93 14 : keys.add("atoms","POSITIONS1","the positions to use for the molecules specified using the first argument");
94 14 : keys.add("atoms","POSITIONS2","the positions to use for the molecules specified using the second argument");
95 28 : keys.setValueDescription("matrix","the matrix of torsions between the two vectors of input directors");
96 14 : }
97 :
98 7 : TorsionsMatrix::TorsionsMatrix(const ActionOptions&ao):
99 : Action(ao),
100 7 : ActionWithMatrix(ao) {
101 7 : if( getNumberOfArguments()!=2 ) {
102 0 : error("should be two arguments to this action, a matrix and a vector");
103 : }
104 7 : if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) {
105 0 : error("first argument to this action should be a matrix");
106 : }
107 7 : if( getPntrToArgument(1)->getRank()!=2 || getPntrToArgument(1)->hasDerivatives() ) {
108 0 : error("second argument to this action should be a matrix");
109 : }
110 7 : if( getPntrToArgument(0)->getShape()[1]!=3 || getPntrToArgument(1)->getShape()[0]!=3 ) {
111 0 : error("number of columns in first matrix and number of rows in second matrix should equal 3");
112 : }
113 :
114 : std::vector<AtomNumber> atoms_a;
115 14 : parseAtomList("POSITIONS1", atoms_a );
116 7 : if( atoms_a.size()!=getPntrToArgument(0)->getShape()[0] ) {
117 0 : error("mismatch between number of atoms specified using POSITIONS1 and number of arguments in vector input");
118 : }
119 7 : log.printf(" using positions of these atoms for vectors in first matrix \n");
120 933 : for(unsigned int i=0; i<atoms_a.size(); ++i) {
121 926 : if ( (i+1) % 25 == 0 ) {
122 36 : log.printf(" \n");
123 : }
124 926 : log.printf(" %d", atoms_a[i].serial());
125 : }
126 7 : log.printf("\n");
127 : std::vector<AtomNumber> atoms_b;
128 14 : parseAtomList("POSITIONS2", atoms_b );
129 7 : if( atoms_b.size()!=getPntrToArgument(1)->getShape()[1] ) {
130 0 : error("mismatch between number of atoms specified using POSITIONS2 and number of arguments in vector input");
131 : }
132 7 : log.printf(" using positions of these atoms for vectors in second matrix \n");
133 1225 : for(unsigned i=0; i<atoms_b.size(); ++i) {
134 1218 : if ( (i+1) % 25 == 0 ) {
135 48 : log.printf(" \n");
136 : }
137 1218 : log.printf(" %d", atoms_b[i].serial());
138 1218 : atoms_a.push_back( atoms_b[i] );
139 : }
140 7 : log.printf("\n");
141 7 : requestAtoms( atoms_a, false );
142 :
143 7 : std::vector<unsigned> shape(2);
144 7 : shape[0]=getPntrToArgument(0)->getShape()[0];
145 7 : shape[1]=getPntrToArgument(1)->getShape()[1];
146 7 : addValue( shape );
147 14 : setPeriodic("-pi","pi");
148 7 : nderivatives = buildArgumentStore(0) + 3*getNumberOfAtoms() + 9;
149 7 : std::string headstr=getFirstActionInChain()->getLabel();
150 7 : stored_matrix1 = getPntrToArgument(0)->ignoreStoredValue( headstr );
151 7 : stored_matrix2 = getPntrToArgument(1)->ignoreStoredValue( headstr );
152 7 : }
153 :
154 10 : unsigned TorsionsMatrix::getNumberOfDerivatives() {
155 10 : return nderivatives;
156 : }
157 :
158 2 : void TorsionsMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
159 2 : unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(1)->getShape()[1];
160 2 : if( indices.size()!=size_v+1 ) {
161 1 : indices.resize( size_v+1 );
162 : }
163 6 : for(unsigned i=0; i<size_v; ++i) {
164 4 : indices[i+1] = start_n + i;
165 : }
166 : myvals.setSplitIndex( size_v + 1 );
167 2 : }
168 :
169 1240496 : void TorsionsMatrix::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
170 1240496 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream(), ind2=index2;
171 1240496 : if( index2>=getPntrToArgument(0)->getShape()[0] ) {
172 26404 : ind2 = index2 - getPntrToArgument(0)->getShape()[0];
173 : }
174 :
175 1240496 : Vector v1, v2, dv1, dv2, dconn;
176 : // Compute the distance connecting the two centers
177 1240496 : Vector conn=pbcDistance( getPosition(index1), getPosition(index2) );
178 1240496 : if( conn.modulo2()<epsilon ) {
179 1239262 : return;
180 : }
181 :
182 : // Get the two vectors
183 4961624 : for(unsigned i=0; i<3; ++i) {
184 3721218 : v1[i] = getElementOfMatrixArgument( 0, index1, i, myvals );
185 3721218 : v2[i] = getElementOfMatrixArgument( 1, i, ind2, myvals );
186 : }
187 : // Evaluate angle
188 : Torsion t;
189 1240406 : double angle = t.compute( v1, conn, v2, dv1, dconn, dv2 );
190 1240406 : myvals.addValue( ostrn, angle );
191 :
192 1240406 : if( doNotCalculateDerivatives() ) {
193 : return;
194 : }
195 :
196 : // Add the derivatives on the matrices
197 4936 : for(unsigned i=0; i<3; ++i) {
198 3702 : addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, i, dv1[i], myvals );
199 3702 : addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, i, ind2, dv2[i], myvals );
200 : }
201 : // And derivatives on positions
202 1234 : unsigned narg_derivatives = getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
203 4936 : for(unsigned i=0; i<3; ++i) {
204 3702 : myvals.addDerivative( ostrn, narg_derivatives + 3*index1+i, -dconn[i] );
205 3702 : myvals.addDerivative( ostrn, narg_derivatives + 3*index2+i, dconn[i] );
206 3702 : myvals.updateIndex( ostrn, narg_derivatives + 3*index1+i );
207 3702 : myvals.updateIndex( ostrn, narg_derivatives + 3*index2+i );
208 : }
209 : //And virial
210 1234 : Tensor vir( -extProduct( conn, dconn ) );
211 1234 : unsigned virbase = narg_derivatives + 3*getNumberOfAtoms();
212 4936 : for(unsigned i=0; i<3; ++i)
213 14808 : for(unsigned j=0; j<3; ++j ) {
214 11106 : myvals.addDerivative( ostrn, virbase+3*i+j, vir(i,j) );
215 11106 : myvals.updateIndex( ostrn, virbase+3*i+j );
216 : }
217 : }
218 :
219 5366 : void TorsionsMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
220 5366 : if( doNotCalculateDerivatives() || !matrixChainContinues() ) {
221 : return ;
222 : }
223 :
224 178 : unsigned mat1s = 3*ival, ss = getPntrToArgument(1)->getShape()[1];
225 178 : unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
226 178 : unsigned narg_derivatives = getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
227 : std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
228 : unsigned ntwo_atoms = myvals.getSplitIndex();
229 712 : for(unsigned j=0; j<3; ++j) {
230 534 : matrix_indices[nmat_ind] = mat1s + j;
231 534 : nmat_ind++;
232 534 : matrix_indices[nmat_ind] = narg_derivatives + mat1s + j;
233 534 : nmat_ind++;
234 4242 : for(unsigned i=1; i<ntwo_atoms; ++i) {
235 3708 : unsigned ind2 = indices[i];
236 3708 : if( ind2>=getPntrToArgument(0)->getShape()[0] ) {
237 12 : ind2 = indices[i] - getPntrToArgument(0)->getShape()[0];
238 : }
239 3708 : matrix_indices[nmat_ind] = arg_deriv_starts[1] + j*ss + ind2;
240 3708 : nmat_ind++;
241 3708 : matrix_indices[nmat_ind] = narg_derivatives + 3*indices[i] + j;
242 3708 : nmat_ind++;
243 : }
244 : }
245 178 : unsigned base = narg_derivatives + 3*getNumberOfAtoms();
246 1780 : for(unsigned j=0; j<9; ++j) {
247 1602 : matrix_indices[nmat_ind] = base + j;
248 1602 : nmat_ind++;
249 : }
250 : myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
251 : }
252 :
253 : }
254 : }
|
