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 : \par Examples
31 :
32 : */
33 : //+ENDPLUMEDOC
34 :
35 : namespace PLMD {
36 : namespace adjmat {
37 :
38 : class TorsionsMatrix : public ActionWithMatrix {
39 : private:
40 : unsigned nderivatives;
41 : bool stored_matrix1, stored_matrix2;
42 : public:
43 : static void registerKeywords( Keywords& keys );
44 : explicit TorsionsMatrix(const ActionOptions&);
45 : unsigned getNumberOfDerivatives();
46 0 : unsigned getNumberOfColumns() const override {
47 0 : return getConstPntrToComponent(0)->getShape()[1];
48 : }
49 : void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
50 : void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
51 : void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
52 : };
53 :
54 : PLUMED_REGISTER_ACTION(TorsionsMatrix,"TORSIONS_MATRIX")
55 :
56 14 : void TorsionsMatrix::registerKeywords( Keywords& keys ) {
57 14 : ActionWithMatrix::registerKeywords(keys);
58 14 : keys.use("ARG");
59 28 : keys.add("atoms","POSITIONS1","the positions to use for the molecules specified using the first argument");
60 28 : keys.add("atoms","POSITIONS2","the positions to use for the molecules specified using the second argument");
61 14 : keys.setValueDescription("the matrix of torsions between the two vectors of input directors");
62 14 : }
63 :
64 7 : TorsionsMatrix::TorsionsMatrix(const ActionOptions&ao):
65 : Action(ao),
66 7 : ActionWithMatrix(ao) {
67 7 : if( getNumberOfArguments()!=2 ) {
68 0 : error("should be two arguments to this action, a matrix and a vector");
69 : }
70 7 : if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) {
71 0 : error("first argument to this action should be a matrix");
72 : }
73 7 : if( getPntrToArgument(1)->getRank()!=2 || getPntrToArgument(1)->hasDerivatives() ) {
74 0 : error("second argument to this action should be a matrix");
75 : }
76 7 : if( getPntrToArgument(0)->getShape()[1]!=3 || getPntrToArgument(1)->getShape()[0]!=3 ) {
77 0 : error("number of columns in first matrix and number of rows in second matrix should equal 3");
78 : }
79 :
80 : std::vector<AtomNumber> atoms_a;
81 14 : parseAtomList("POSITIONS1", atoms_a );
82 7 : if( atoms_a.size()!=getPntrToArgument(0)->getShape()[0] ) {
83 0 : error("mismatch between number of atoms specified using POSITIONS1 and number of arguments in vector input");
84 : }
85 7 : log.printf(" using positions of these atoms for vectors in first matrix \n");
86 933 : for(unsigned int i=0; i<atoms_a.size(); ++i) {
87 926 : if ( (i+1) % 25 == 0 ) {
88 36 : log.printf(" \n");
89 : }
90 926 : log.printf(" %d", atoms_a[i].serial());
91 : }
92 7 : log.printf("\n");
93 : std::vector<AtomNumber> atoms_b;
94 14 : parseAtomList("POSITIONS2", atoms_b );
95 7 : if( atoms_b.size()!=getPntrToArgument(1)->getShape()[1] ) {
96 0 : error("mismatch between number of atoms specified using POSITIONS2 and number of arguments in vector input");
97 : }
98 7 : log.printf(" using positions of these atoms for vectors in second matrix \n");
99 1225 : for(unsigned i=0; i<atoms_b.size(); ++i) {
100 1218 : if ( (i+1) % 25 == 0 ) {
101 48 : log.printf(" \n");
102 : }
103 1218 : log.printf(" %d", atoms_b[i].serial());
104 1218 : atoms_a.push_back( atoms_b[i] );
105 : }
106 7 : log.printf("\n");
107 7 : requestAtoms( atoms_a, false );
108 :
109 7 : std::vector<unsigned> shape(2);
110 7 : shape[0]=getPntrToArgument(0)->getShape()[0];
111 7 : shape[1]=getPntrToArgument(1)->getShape()[1];
112 7 : addValue( shape );
113 14 : setPeriodic("-pi","pi");
114 7 : nderivatives = buildArgumentStore(0) + 3*getNumberOfAtoms() + 9;
115 7 : std::string headstr=getFirstActionInChain()->getLabel();
116 7 : stored_matrix1 = getPntrToArgument(0)->ignoreStoredValue( headstr );
117 7 : stored_matrix2 = getPntrToArgument(1)->ignoreStoredValue( headstr );
118 7 : }
119 :
120 10 : unsigned TorsionsMatrix::getNumberOfDerivatives() {
121 10 : return nderivatives;
122 : }
123 :
124 2 : void TorsionsMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
125 2 : unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(1)->getShape()[1];
126 2 : if( indices.size()!=size_v+1 ) {
127 1 : indices.resize( size_v+1 );
128 : }
129 6 : for(unsigned i=0; i<size_v; ++i) {
130 4 : indices[i+1] = start_n + i;
131 : }
132 : myvals.setSplitIndex( size_v + 1 );
133 2 : }
134 :
135 1240496 : void TorsionsMatrix::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
136 1240496 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream(), ind2=index2;
137 1240496 : if( index2>=getPntrToArgument(0)->getShape()[0] ) {
138 26404 : ind2 = index2 - getPntrToArgument(0)->getShape()[0];
139 : }
140 :
141 1240496 : Vector v1, v2, dv1, dv2, dconn;
142 : // Compute the distance connecting the two centers
143 1240496 : Vector conn=pbcDistance( getPosition(index1), getPosition(index2) );
144 1240496 : if( conn.modulo2()<epsilon ) {
145 1239262 : return;
146 : }
147 :
148 : // Get the two vectors
149 4961624 : for(unsigned i=0; i<3; ++i) {
150 3721218 : v1[i] = getElementOfMatrixArgument( 0, index1, i, myvals );
151 3721218 : v2[i] = getElementOfMatrixArgument( 1, i, ind2, myvals );
152 : }
153 : // Evaluate angle
154 : Torsion t;
155 1240406 : double angle = t.compute( v1, conn, v2, dv1, dconn, dv2 );
156 1240406 : myvals.addValue( ostrn, angle );
157 :
158 1240406 : if( doNotCalculateDerivatives() ) {
159 : return;
160 : }
161 :
162 : // Add the derivatives on the matrices
163 4936 : for(unsigned i=0; i<3; ++i) {
164 3702 : addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, i, dv1[i], myvals );
165 3702 : addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, i, ind2, dv2[i], myvals );
166 : }
167 : // And derivatives on positions
168 1234 : unsigned narg_derivatives = getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
169 4936 : for(unsigned i=0; i<3; ++i) {
170 3702 : myvals.addDerivative( ostrn, narg_derivatives + 3*index1+i, -dconn[i] );
171 3702 : myvals.addDerivative( ostrn, narg_derivatives + 3*index2+i, dconn[i] );
172 3702 : myvals.updateIndex( ostrn, narg_derivatives + 3*index1+i );
173 3702 : myvals.updateIndex( ostrn, narg_derivatives + 3*index2+i );
174 : }
175 : //And virial
176 1234 : Tensor vir( -extProduct( conn, dconn ) );
177 1234 : unsigned virbase = narg_derivatives + 3*getNumberOfAtoms();
178 4936 : for(unsigned i=0; i<3; ++i)
179 14808 : for(unsigned j=0; j<3; ++j ) {
180 11106 : myvals.addDerivative( ostrn, virbase+3*i+j, vir(i,j) );
181 11106 : myvals.updateIndex( ostrn, virbase+3*i+j );
182 : }
183 : }
184 :
185 5366 : void TorsionsMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
186 5366 : if( doNotCalculateDerivatives() || !matrixChainContinues() ) {
187 : return ;
188 : }
189 :
190 178 : unsigned mat1s = 3*ival, ss = getPntrToArgument(1)->getShape()[1];
191 178 : unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
192 178 : unsigned narg_derivatives = getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
193 : std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
194 : unsigned ntwo_atoms = myvals.getSplitIndex();
195 712 : for(unsigned j=0; j<3; ++j) {
196 534 : matrix_indices[nmat_ind] = mat1s + j;
197 534 : nmat_ind++;
198 534 : matrix_indices[nmat_ind] = narg_derivatives + mat1s + j;
199 534 : nmat_ind++;
200 4242 : for(unsigned i=1; i<ntwo_atoms; ++i) {
201 3708 : unsigned ind2 = indices[i];
202 3708 : if( ind2>=getPntrToArgument(0)->getShape()[0] ) {
203 12 : ind2 = indices[i] - getPntrToArgument(0)->getShape()[0];
204 : }
205 3708 : matrix_indices[nmat_ind] = arg_deriv_starts[1] + j*ss + ind2;
206 3708 : nmat_ind++;
207 3708 : matrix_indices[nmat_ind] = narg_derivatives + 3*indices[i] + j;
208 3708 : nmat_ind++;
209 : }
210 : }
211 178 : unsigned base = narg_derivatives + 3*getNumberOfAtoms();
212 1780 : for(unsigned j=0; j<9; ++j) {
213 1602 : matrix_indices[nmat_ind] = base + j;
214 1602 : nmat_ind++;
215 : }
216 : myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
217 : }
218 :
219 : }
220 : }
|
