Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2018 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 ?? so two atoms are adjacent if they are within a cutoff:
48 :
49 : \verbatim
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 : \endverbatim
55 :
56 : This example input calculates the 14 SPRINT coordinates foa 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 : \verbatim
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 : \endverbatim
75 :
76 : */
77 : //+ENDPLUMEDOC
78 :
79 :
80 : namespace PLMD {
81 : namespace adjmat {
82 :
83 2 : 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();
104 : /// Sprint needs its only apply routine as it creates values
105 : void apply();
106 : };
107 :
108 2524 : PLUMED_REGISTER_ACTION(Sprint,"SPRINT")
109 :
110 2 : void Sprint::registerKeywords( Keywords& keys ) {
111 2 : ActionWithInputMatrix::registerKeywords( keys );
112 2 : componentsAreNotOptional(keys);
113 : 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 labelled <em>label</em>.coord-1, "
115 2 : "the second smallest will be labelleled <em>label</em>.coord-1 and so on");
116 2 : }
117 :
118 1 : Sprint::Sprint(const ActionOptions&ao):
119 : Action(ao),
120 : ActionWithInputMatrix(ao),
121 1 : eigvals( getNumberOfNodes() ),
122 1 : maxeig( getNumberOfNodes() ),
123 : thematrix( getNumberOfNodes(), getNumberOfNodes() ),
124 3 : eigenvecs( getNumberOfNodes(), getNumberOfNodes() )
125 : {
126 : // Check on setup
127 : // if( getNumberOfVessels()!=1 ) error("there should be no vessel keywords");
128 : // Check for bad colvar input ( we are going to get rid of this because we are going to have input adjacency matrix in future )
129 : // for(unsigned i=0;i<getNumberOfAtomGroups();++i){
130 : // /// Check me GAT
131 : // // if( !getBaseMultiColvar(i)->hasDifferentiableOrientation() ) error("cannot use multicolvar of type " + getBaseMultiColvar(i)->getName() );
132 : // }
133 :
134 1 : if( !getAdjacencyVessel()->isSymmetric() ) error("input contact matrix is not symmetric");
135 1 : std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
136 :
137 : // Create all the values
138 1 : sqrtn = sqrt( static_cast<double>( getNumberOfNodes() ) );
139 15 : for(unsigned i=0; i<getNumberOfNodes(); ++i) {
140 14 : std::string num; Tools::convert(i,num);
141 14 : addComponentWithDerivatives("coord-"+num);
142 14 : componentIsNotPeriodic("coord-"+num);
143 14 : getPntrToComponent(i)->resizeDerivatives( getNumberOfDerivatives() );
144 14 : }
145 :
146 : // Setup the dynamic list to hold all the tasks
147 1 : unsigned ntriangle = 0.5*getNumberOfNodes()*(getNumberOfNodes()-1);
148 1 : for(unsigned i=0; i<ntriangle; ++i) active_elements.addIndexToList( i );
149 1 : }
150 :
151 17 : void Sprint::calculate() {
152 : // Get the adjacency matrix
153 17 : getAdjacencyVessel()->retrieveMatrix( active_elements, thematrix );
154 : // Diagonalize it
155 17 : diagMat( thematrix, eigvals, eigenvecs );
156 : // Get the maximum eigevalue
157 17 : double lambda = eigvals[ getNumberOfNodes()-1 ];
158 : // Get the corresponding eigenvector
159 255 : for(unsigned j=0; j<maxeig.size(); ++j) {
160 238 : maxeig[j].first = fabs( eigenvecs( getNumberOfNodes()-1, j ) );
161 238 : maxeig[j].second = j;
162 : // Must make all components of principle eigenvector +ve
163 238 : eigenvecs( getNumberOfNodes()-1, j ) = maxeig[j].first;
164 : }
165 :
166 : // Reorder each block of eigevectors
167 17 : unsigned startnum=0;
168 51 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
169 34 : unsigned nthis = getNumberOfAtomsInGroup(j);
170 : // Sort into ascending order
171 34 : std::sort( maxeig.begin() + startnum, maxeig.begin() + startnum + nthis );
172 : // Used so we can do sorting in blocks
173 34 : startnum += nthis;
174 : }
175 : // Set the sprint coordinates
176 255 : for(unsigned icomp=0; icomp<getNumberOfComponents(); ++icomp) {
177 238 : getPntrToComponent(icomp)->set( sqrtn*lambda*maxeig[icomp].first );
178 : }
179 :
180 : // Parallelism
181 : unsigned rank, stride;
182 17 : if( serialCalculation() ) { stride=1; rank=0; }
183 17 : else { rank=comm.Get_rank(); stride=comm.Get_size(); }
184 :
185 : // Derivatives
186 17 : MultiValue myvals( 2, getNumberOfDerivatives() );
187 34 : Matrix<double> mymat_ders( getNumberOfComponents(), getNumberOfDerivatives() );
188 : // std::vector<unsigned> catoms(2);
189 17 : unsigned nval = getNumberOfNodes(); mymat_ders=0;
190 1564 : for(unsigned i=rank; i<active_elements.getNumberActive(); i+=stride) {
191 1547 : unsigned j, k; getAdjacencyVessel()->getMatrixIndices( active_elements[i], j, k );
192 1547 : double tmp1 = 2 * eigenvecs(nval-1,j)*eigenvecs(nval-1,k);
193 23205 : for(unsigned icomp=0; icomp<getNumberOfComponents(); ++icomp) {
194 21658 : double tmp2 = 0.;
195 303212 : for(unsigned n=0; n<nval-1; ++n) { // Need care on following line
196 281554 : tmp2 += eigenvecs(n,maxeig[icomp].second) * ( eigenvecs(n,j)*eigenvecs(nval-1,k) + eigenvecs(n,k)*eigenvecs(nval-1,j) ) / ( lambda - eigvals[n] );
197 : }
198 21658 : double prefactor=sqrtn*( tmp1*maxeig[icomp].first + tmp2*lambda );
199 21658 : getAdjacencyVessel()->retrieveDerivatives( active_elements[i], false, myvals );
200 346528 : for(unsigned jd=0; jd<myvals.getNumberActive(); ++jd) {
201 324870 : unsigned ider=myvals.getActiveIndex(jd);
202 324870 : mymat_ders( icomp, ider ) += prefactor*myvals.getDerivative( 1, ider );
203 : }
204 : }
205 : }
206 17 : if( !serialCalculation() ) comm.Sum( mymat_ders );
207 :
208 255 : for(unsigned j=0; j<getNumberOfComponents(); ++j) {
209 238 : Value* val=getPntrToComponent(j);
210 238 : for(unsigned i=0; i<getNumberOfDerivatives(); ++i) val->addDerivative( i, mymat_ders(j,i) );
211 17 : }
212 17 : }
213 :
214 17 : void Sprint::apply() {
215 17 : std::vector<Vector>& f(modifyForces());
216 17 : Tensor& v(modifyVirial());
217 17 : unsigned nat=getNumberOfAtoms();
218 :
219 17 : std::vector<double> forces( 3*getNumberOfAtoms() + 9 );
220 255 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
221 238 : if( getPntrToComponent(i)->applyForce( forces ) ) {
222 0 : for(unsigned j=0; j<nat; ++j) {
223 0 : f[j][0]+=forces[3*j+0];
224 0 : f[j][1]+=forces[3*j+1];
225 0 : f[j][2]+=forces[3*j+2];
226 : }
227 0 : v(0,0)+=forces[3*nat+0];
228 0 : v(0,1)+=forces[3*nat+1];
229 0 : v(0,2)+=forces[3*nat+2];
230 0 : v(1,0)+=forces[3*nat+3];
231 0 : v(1,1)+=forces[3*nat+4];
232 0 : v(1,2)+=forces[3*nat+5];
233 0 : v(2,0)+=forces[3*nat+6];
234 0 : v(2,1)+=forces[3*nat+7];
235 0 : v(2,2)+=forces[3*nat+8];
236 : }
237 17 : }
238 17 : }
239 :
240 : }
241 2523 : }
|
