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 { return getConstPntrToComponent(0)->getShape()[1]; }
47 : void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
48 : void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
49 : void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
50 : };
51 :
52 : PLUMED_REGISTER_ACTION(TorsionsMatrix,"TORSIONS_MATRIX")
53 :
54 14 : void TorsionsMatrix::registerKeywords( Keywords& keys ) {
55 14 : ActionWithMatrix::registerKeywords(keys); keys.use("ARG");
56 28 : keys.add("atoms","POSITIONS1","the positions to use for the molecules specified using the first argument");
57 28 : keys.add("atoms","POSITIONS2","the positions to use for the molecules specified using the second argument");
58 14 : keys.setValueDescription("the matrix of torsions between the two vectors of input directors");
59 14 : }
60 :
61 7 : TorsionsMatrix::TorsionsMatrix(const ActionOptions&ao):
62 : Action(ao),
63 7 : ActionWithMatrix(ao)
64 : {
65 7 : if( getNumberOfArguments()!=2 ) error("should be two arguments to this action, a matrix and a vector");
66 7 : if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) error("first argument to this action should be a matrix");
67 7 : if( getPntrToArgument(1)->getRank()!=2 || getPntrToArgument(1)->hasDerivatives() ) error("second argument to this action should be a matrix");
68 7 : if( getPntrToArgument(0)->getShape()[1]!=3 || getPntrToArgument(1)->getShape()[0]!=3 ) error("number of columns in first matrix and number of rows in second matrix should equal 3");
69 :
70 14 : std::vector<AtomNumber> atoms_a; parseAtomList("POSITIONS1", atoms_a );
71 7 : if( atoms_a.size()!=getPntrToArgument(0)->getShape()[0] ) error("mismatch between number of atoms specified using POSITIONS1 and number of arguments in vector input");
72 7 : log.printf(" using positions of these atoms for vectors in first matrix \n");
73 933 : for(unsigned int i=0; i<atoms_a.size(); ++i) {
74 926 : if ( (i+1) % 25 == 0 ) log.printf(" \n");
75 926 : log.printf(" %d", atoms_a[i].serial());
76 : }
77 14 : log.printf("\n"); std::vector<AtomNumber> atoms_b; parseAtomList("POSITIONS2", atoms_b );
78 7 : if( atoms_b.size()!=getPntrToArgument(1)->getShape()[1] ) error("mismatch between number of atoms specified using POSITIONS2 and number of arguments in vector input");
79 7 : log.printf(" using positions of these atoms for vectors in second matrix \n");
80 1225 : for(unsigned i=0; i<atoms_b.size(); ++i) {
81 1218 : if ( (i+1) % 25 == 0 ) log.printf(" \n");
82 1218 : log.printf(" %d", atoms_b[i].serial()); atoms_a.push_back( atoms_b[i] );
83 : }
84 7 : log.printf("\n"); requestAtoms( atoms_a, false );
85 :
86 7 : std::vector<unsigned> shape(2); shape[0]=getPntrToArgument(0)->getShape()[0]; shape[1]=getPntrToArgument(1)->getShape()[1];
87 14 : addValue( shape ); setPeriodic("-pi","pi"); nderivatives = buildArgumentStore(0) + 3*getNumberOfAtoms() + 9;
88 7 : std::string headstr=getFirstActionInChain()->getLabel();
89 7 : stored_matrix1 = getPntrToArgument(0)->ignoreStoredValue( headstr );
90 7 : stored_matrix2 = getPntrToArgument(1)->ignoreStoredValue( headstr );
91 7 : }
92 :
93 10 : unsigned TorsionsMatrix::getNumberOfDerivatives() {
94 10 : return nderivatives;
95 : }
96 :
97 2 : void TorsionsMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
98 2 : unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(1)->getShape()[1];
99 2 : if( indices.size()!=size_v+1 ) indices.resize( size_v+1 );
100 6 : for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + i;
101 : myvals.setSplitIndex( size_v + 1 );
102 2 : }
103 :
104 1240496 : void TorsionsMatrix::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
105 1240496 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream(), ind2=index2;
106 1240496 : if( index2>=getPntrToArgument(0)->getShape()[0] ) ind2 = index2 - getPntrToArgument(0)->getShape()[0];
107 :
108 1240496 : Vector v1, v2, dv1, dv2, dconn;
109 : // Compute the distance connecting the two centers
110 1240496 : Vector conn=pbcDistance( getPosition(index1), getPosition(index2) );
111 2479668 : if( conn.modulo2()<epsilon ) return;
112 :
113 : // Get the two vectors
114 4961624 : for(unsigned i=0; i<3; ++i) {
115 3721218 : v1[i] = getElementOfMatrixArgument( 0, index1, i, myvals );
116 3721218 : v2[i] = getElementOfMatrixArgument( 1, i, ind2, myvals );
117 : }
118 : // Evaluate angle
119 1240406 : Torsion t; double angle = t.compute( v1, conn, v2, dv1, dconn, dv2 );
120 1240406 : myvals.addValue( ostrn, angle );
121 :
122 1240406 : if( doNotCalculateDerivatives() ) return;
123 :
124 : // Add the derivatives on the matrices
125 4936 : for(unsigned i=0; i<3; ++i) {
126 3702 : addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, i, dv1[i], myvals );
127 3702 : addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, i, ind2, dv2[i], myvals );
128 : }
129 : // And derivatives on positions
130 1234 : unsigned narg_derivatives = getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
131 4936 : for(unsigned i=0; i<3; ++i) {
132 3702 : myvals.addDerivative( ostrn, narg_derivatives + 3*index1+i, -dconn[i] ); myvals.addDerivative( ostrn, narg_derivatives + 3*index2+i, dconn[i] );
133 3702 : myvals.updateIndex( ostrn, narg_derivatives + 3*index1+i ); myvals.updateIndex( ostrn, narg_derivatives + 3*index2+i );
134 : }
135 : //And virial
136 1234 : Tensor vir( -extProduct( conn, dconn ) ); unsigned virbase = narg_derivatives + 3*getNumberOfAtoms();
137 16042 : for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j ) { myvals.addDerivative( ostrn, virbase+3*i+j, vir(i,j) ); myvals.updateIndex( ostrn, virbase+3*i+j ); }
138 : }
139 :
140 5366 : void TorsionsMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
141 5366 : if( doNotCalculateDerivatives() || !matrixChainContinues() ) return ;
142 :
143 178 : unsigned mat1s = 3*ival, ss = getPntrToArgument(1)->getShape()[1];
144 178 : unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
145 178 : unsigned narg_derivatives = getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
146 : std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); unsigned ntwo_atoms = myvals.getSplitIndex();
147 712 : for(unsigned j=0; j<3; ++j) {
148 534 : matrix_indices[nmat_ind] = mat1s + j; nmat_ind++;
149 534 : matrix_indices[nmat_ind] = narg_derivatives + mat1s + j; nmat_ind++;
150 4242 : for(unsigned i=1; i<ntwo_atoms; ++i) {
151 3708 : unsigned ind2 = indices[i]; if( ind2>=getPntrToArgument(0)->getShape()[0] ) ind2 = indices[i] - getPntrToArgument(0)->getShape()[0];
152 3708 : matrix_indices[nmat_ind] = arg_deriv_starts[1] + j*ss + ind2; nmat_ind++;
153 3708 : matrix_indices[nmat_ind] = narg_derivatives + 3*indices[i] + j; nmat_ind++;
154 : }
155 : }
156 1780 : unsigned base = narg_derivatives + 3*getNumberOfAtoms(); for(unsigned j=0; j<9; ++j) { matrix_indices[nmat_ind] = base + j; nmat_ind++; }
157 : myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
158 : }
159 :
160 : }
161 : }
|