Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "ActionWithInputMatrix.h"
23 : #include "AdjacencyMatrixVessel.h"
24 : #include "core/ActionRegister.h"
25 :
26 : //+PLUMEDOC MATRIXF SPRINT
27 : /*
28 : Calculate SPRINT topological variables from an adjacency matrix.
29 :
30 : The SPRINT topological variables are calculated from the largest eigenvalue, \f$\lambda\f$ of
31 : an \f$n\times n\f$ adjacency matrix and its corresponding eigenvector, \f$\mathbf{V}\f$, using:
32 :
33 : \f[
34 : s_i = \sqrt{n} \lambda v_i
35 : \f]
36 :
37 : You can use different quantities to measure whether or not two given atoms/molecules are
38 : adjacent or not in the adjacency matrix. The simplest measure of adjacency is is whether
39 : two atoms/molecules are within some cutoff of each other. Further complexity can be added by
40 : insisting that two molecules are adjacent if they are within a certain distance of each
41 : other and if they have similar orientations.
42 :
43 : \par Examples
44 :
45 : This example input calculates the 7 SPRINT coordinates for a 7 atom cluster of Lennard-Jones
46 : atoms and prints their values to a file. In this input the SPRINT coordinates are calculated
47 : in the manner described in \cite Pietrucci2011 so two atoms are adjacent if they are within a cutoff:
48 :
49 : \plumedfile
50 : DENSITY SPECIES=1-7 LABEL=d1
51 : CONTACT_MATRIX ATOMS=d1 SWITCH={RATIONAL R_0=0.1} LABEL=mat
52 : SPRINT MATRIX=mat LABEL=ss
53 : PRINT ARG=ss.* FILE=colvar
54 : \endplumedfile
55 :
56 : This example input calculates the 14 SPRINT coordinates for a molecule composed of 7 hydrogen and
57 : 7 carbon atoms. Once again two atoms are adjacent if they are within a cutoff:
58 :
59 : \plumedfile
60 : DENSITY SPECIES=1-7 LABEL=c
61 : DENSITY SPECIES=8-14 LABEL=h
62 :
63 : CONTACT_MATRIX ...
64 : ATOMS=c,h
65 : SWITCH11={RATIONAL R_0=2.6 NN=6 MM=12}
66 : SWITCH12={RATIONAL R_0=2.2 NN=6 MM=12}
67 : SWITCH22={RATIONAL R_0=2.2 NN=6 MM=12}
68 : LABEL=mat
69 : ... CONTACT_MATRIX
70 :
71 : SPRINT MATRIX=mat LABEL=ss
72 :
73 : PRINT ARG=ss.* FILE=colvar
74 : \endplumedfile
75 :
76 : */
77 : //+ENDPLUMEDOC
78 :
79 :
80 : namespace PLMD {
81 : namespace adjmat {
82 :
83 : class Sprint : public ActionWithInputMatrix {
84 : private:
85 : /// Square root of number of atoms
86 : double sqrtn;
87 : /// Vector that stores eigenvalues
88 : std::vector<double> eigvals;
89 : /// This is used to speed up the calculation of derivatives
90 : DynamicList<unsigned> active_elements;
91 : /// Vector that stores max eigenvector
92 : std::vector< std::pair<double,int> > maxeig;
93 : /// Adjacency matrix
94 : Matrix<double> thematrix;
95 : /// Matrix that stores eigenvectors
96 : Matrix<double> eigenvecs;
97 : public:
98 : /// Create manual
99 : static void registerKeywords( Keywords& keys );
100 : /// Constructor
101 : explicit Sprint(const ActionOptions&);
102 : /// Do the matrix calculation
103 : void calculate() override;
104 : /// Sprint needs its only apply routine as it creates values
105 : void apply() override;
106 : };
107 :
108 13787 : PLUMED_REGISTER_ACTION(Sprint,"SPRINT")
109 :
110 5 : void Sprint::registerKeywords( Keywords& keys ) {
111 5 : ActionWithInputMatrix::registerKeywords( keys );
112 5 : componentsAreNotOptional(keys);
113 10 : keys.addOutputComponent("coord","default","all \\f$n\\f$ sprint coordinates are calculated and then stored in increasing order. "
114 : "the smallest sprint coordinate will be labeled <em>label</em>.coord-0, "
115 : "the second smallest will be labelled <em>label</em>.coord-1 and so on");
116 5 : }
117 :
118 1 : Sprint::Sprint(const ActionOptions&ao):
119 : Action(ao),
120 : ActionWithInputMatrix(ao),
121 1 : eigvals( getNumberOfNodes() ),
122 1 : maxeig( getNumberOfNodes() ),
123 1 : thematrix( getNumberOfNodes(), getNumberOfNodes() ),
124 3 : eigenvecs( getNumberOfNodes(), getNumberOfNodes() ) {
125 : // Check on setup
126 : // if( getNumberOfVessels()!=1 ) error("there should be no vessel keywords");
127 : // Check for bad colvar input ( we are going to get rid of this because we are going to have input adjacency matrix in future )
128 : // for(unsigned i=0;i<getNumberOfAtomGroups();++i){
129 : // /// Check me GAT
130 : // // if( !getBaseMultiColvar(i)->hasDifferentiableOrientation() ) error("cannot use multicolvar of type " + getBaseMultiColvar(i)->getName() );
131 : // }
132 :
133 1 : if( !getAdjacencyVessel()->isSymmetric() ) {
134 0 : error("input contact matrix is not symmetric");
135 : }
136 : std::vector<AtomNumber> fake_atoms;
137 1 : setupMultiColvarBase( fake_atoms );
138 :
139 : // Create all the values
140 1 : sqrtn = std::sqrt( static_cast<double>( getNumberOfNodes() ) );
141 15 : for(unsigned i=0; i<getNumberOfNodes(); ++i) {
142 : std::string num;
143 14 : Tools::convert(i,num);
144 14 : addComponentWithDerivatives("coord-"+num);
145 14 : componentIsNotPeriodic("coord-"+num);
146 14 : getPntrToComponent(i)->resizeDerivatives( getNumberOfDerivatives() );
147 : }
148 :
149 : // Setup the dynamic list to hold all the tasks
150 1 : unsigned ntriangle = 0.5*getNumberOfNodes()*(getNumberOfNodes()-1);
151 92 : for(unsigned i=0; i<ntriangle; ++i) {
152 91 : active_elements.addIndexToList( i );
153 : }
154 1 : }
155 :
156 17 : void Sprint::calculate() {
157 : // Get the adjacency matrix
158 17 : getAdjacencyVessel()->retrieveMatrix( active_elements, thematrix );
159 : // Diagonalize it
160 17 : diagMat( thematrix, eigvals, eigenvecs );
161 : // Get the maximum eigevalue
162 17 : double lambda = eigvals[ getNumberOfNodes()-1 ];
163 : // Get the corresponding eigenvector
164 255 : for(unsigned j=0; j<maxeig.size(); ++j) {
165 238 : maxeig[j].first = std::fabs( eigenvecs( getNumberOfNodes()-1, j ) );
166 238 : maxeig[j].second = j;
167 : // Must make all components of principle eigenvector +ve
168 238 : eigenvecs( getNumberOfNodes()-1, j ) = maxeig[j].first;
169 : }
170 :
171 : // Reorder each block of eigevectors
172 : unsigned startnum=0;
173 51 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
174 34 : unsigned nthis = getNumberOfAtomsInGroup(j);
175 : // Sort into ascending order
176 34 : std::sort( maxeig.begin() + startnum, maxeig.begin() + startnum + nthis );
177 : // Used so we can do sorting in blocks
178 34 : startnum += nthis;
179 : }
180 : // Set the sprint coordinates
181 255 : for(int icomp=0; icomp<getNumberOfComponents(); ++icomp) {
182 238 : getPntrToComponent(icomp)->set( sqrtn*lambda*maxeig[icomp].first );
183 : }
184 :
185 : // Parallelism
186 : unsigned rank, stride;
187 17 : if( serialCalculation() ) {
188 : stride=1;
189 : rank=0;
190 : } else {
191 17 : rank=comm.Get_rank();
192 17 : stride=comm.Get_size();
193 : }
194 :
195 : // Derivatives
196 17 : MultiValue myvals( 2, getNumberOfDerivatives() );
197 17 : Matrix<double> mymat_ders( getNumberOfComponents(), getNumberOfDerivatives() );
198 : // std::vector<unsigned> catoms(2);
199 17 : unsigned nval = getNumberOfNodes();
200 : mymat_ders=0;
201 1564 : for(unsigned i=rank; i<active_elements.getNumberActive(); i+=stride) {
202 : unsigned j, k;
203 1547 : getAdjacencyVessel()->getMatrixIndices( active_elements[i], j, k );
204 1547 : double tmp1 = 2 * eigenvecs(nval-1,j)*eigenvecs(nval-1,k);
205 23205 : for(int icomp=0; icomp<getNumberOfComponents(); ++icomp) {
206 : double tmp2 = 0.;
207 303212 : for(unsigned n=0; n<nval-1; ++n) { // Need care on following line
208 281554 : tmp2 += eigenvecs(n,maxeig[icomp].second) * ( eigenvecs(n,j)*eigenvecs(nval-1,k) + eigenvecs(n,k)*eigenvecs(nval-1,j) ) / ( lambda - eigvals[n] );
209 : }
210 21658 : double prefactor=sqrtn*( tmp1*maxeig[icomp].first + tmp2*lambda );
211 21658 : getAdjacencyVessel()->retrieveDerivatives( active_elements[i], false, myvals );
212 346528 : for(unsigned jd=0; jd<myvals.getNumberActive(); ++jd) {
213 : unsigned ider=myvals.getActiveIndex(jd);
214 324870 : mymat_ders( icomp, ider ) += prefactor*myvals.getDerivative( 1, ider );
215 : }
216 : }
217 : }
218 17 : if( !serialCalculation() ) {
219 17 : comm.Sum( mymat_ders );
220 : }
221 :
222 255 : for(int j=0; j<getNumberOfComponents(); ++j) {
223 238 : Value* val=getPntrToComponent(j);
224 12376 : for(unsigned i=0; i<getNumberOfDerivatives(); ++i) {
225 12138 : val->addDerivative( i, mymat_ders(j,i) );
226 : }
227 : }
228 17 : }
229 :
230 17 : void Sprint::apply() {
231 : std::vector<Vector>& f(modifyForces());
232 : Tensor& v(modifyVirial());
233 : unsigned nat=getNumberOfAtoms();
234 :
235 17 : std::vector<double> forces( 3*getNumberOfAtoms() + 9 );
236 255 : for(int i=0; i<getNumberOfComponents(); ++i) {
237 238 : if( getPntrToComponent(i)->applyForce( forces ) ) {
238 0 : for(unsigned j=0; j<nat; ++j) {
239 0 : f[j][0]+=forces[3*j+0];
240 0 : f[j][1]+=forces[3*j+1];
241 0 : f[j][2]+=forces[3*j+2];
242 : }
243 0 : v(0,0)+=forces[3*nat+0];
244 0 : v(0,1)+=forces[3*nat+1];
245 0 : v(0,2)+=forces[3*nat+2];
246 0 : v(1,0)+=forces[3*nat+3];
247 0 : v(1,1)+=forces[3*nat+4];
248 0 : v(1,2)+=forces[3*nat+5];
249 0 : v(2,0)+=forces[3*nat+6];
250 0 : v(2,1)+=forces[3*nat+7];
251 0 : v(2,2)+=forces[3*nat+8];
252 : }
253 : }
254 17 : }
255 :
256 : }
257 : }
|
