Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2019 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 "AdjacencyMatrixBase.h"
23 : #include "multicolvar/AtomValuePack.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/SwitchingFunction.h"
26 : #include "tools/HistogramBead.h"
27 : #include "tools/Matrix.h"
28 :
29 : //+PLUMEDOC MATRIX TOPOLOGY_MATRIX
30 : /*
31 : Adjacency matrix in which two atoms are adjacent if they are connected topologically
32 :
33 : \par Examples
34 :
35 :
36 : */
37 : //+ENDPLUMEDOC
38 :
39 : namespace PLMD {
40 : namespace adjmat {
41 :
42 18 : class TopologyMatrix : public AdjacencyMatrixBase {
43 : private:
44 : /// The width to use for the kernel density estimation and the
45 : /// sizes of the bins to be used in kernel density estimation
46 : double sigma;
47 : std::string kerneltype;
48 : /// The maximum number of bins that will be used
49 : /// This is calculated based on the dmax of the switching functions
50 : unsigned maxbins;
51 : /// The volume of the cells
52 : Matrix<double> cell_volume;
53 : /// switching function
54 : Matrix<SwitchingFunction> switchingFunction;
55 : Matrix<SwitchingFunction> cylinder_sw;
56 : Matrix<SwitchingFunction> low_sf;
57 : double beadrad, lsfmax;
58 : Matrix<double> binw_mat;
59 : SwitchingFunction threshold_switch;
60 : public:
61 : /// Create manual
62 : static void registerKeywords( Keywords& keys );
63 : /// Constructor
64 : explicit TopologyMatrix(const ActionOptions&);
65 : /// Get the number of quantities that we must compute
66 : unsigned getNumberOfQuantities() const ;
67 : /// Create the ith, ith switching function
68 : void setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc );
69 : /// This actually calculates the value of the contact function
70 : double calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const ;
71 : /// This does nothing
72 : double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
73 : /// Calculate the contribution from one of the atoms in the third element of the pack
74 : void calculateForThreeAtoms( const unsigned& iat, const Vector& d1, const double& d1_len,
75 : HistogramBead& bead, multicolvar::AtomValuePack& myatoms ) const ;
76 : };
77 :
78 6458 : PLUMED_REGISTER_ACTION(TopologyMatrix,"TOPOLOGY_MATRIX")
79 :
80 7 : void TopologyMatrix::registerKeywords( Keywords& keys ) {
81 7 : AdjacencyMatrixBase::registerKeywords( keys );
82 28 : keys.add("atoms","NODES","The list of atoms for which you would like to calculate the contact matrix. The atoms involved must be specified "
83 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
84 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
85 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
86 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
87 28 : keys.add("atoms","ATOMS","");
88 28 : keys.add("numbered","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
89 : "The following provides information on the \\ref switchingfunction that are available. "
90 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
91 28 : keys.add("numbered","RADIUS","");
92 28 : keys.add("numbered","CYLINDER_SWITCH","a switching function on ( r_ij . r_ik - 1 )/r_ij");
93 28 : keys.add("numbered","BIN_SIZE","the size to use for the bins");
94 28 : keys.add("compulsory","DENSITY_THRESHOLD","");
95 28 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
96 35 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
97 28 : keys.add("hidden","FAKE","");
98 7 : }
99 :
100 6 : TopologyMatrix::TopologyMatrix( const ActionOptions& ao ):
101 : Action(ao),
102 : AdjacencyMatrixBase(ao),
103 18 : maxbins(0)
104 : {
105 30 : readMaxThreeSpeciesMatrix("NODES", "FAKE", "FAKE", "ATOMS", true );
106 6 : unsigned nrows, ncols, ndonor_types; retrieveTypeDimensions( nrows, ncols, ndonor_types );
107 18 : switchingFunction.resize( nrows, ncols ); parseConnectionDescriptions("SWITCH",false,ndonor_types);
108 18 : cylinder_sw.resize( nrows, ncols ); parseConnectionDescriptions("RADIUS",false,ndonor_types);
109 18 : low_sf.resize( nrows, ncols ); parseConnectionDescriptions("CYLINDER_SWITCH",false,ndonor_types);
110 12 : binw_mat.resize( nrows, ncols ); cell_volume.resize( nrows, ncols );
111 12 : parseConnectionDescriptions("BIN_SIZE",false,ndonor_types);
112 : // Read in stuff for grid
113 18 : parse("SIGMA",sigma); parse("KERNEL",kerneltype);
114 : // Read in threshold for density cutoff
115 12 : std::string errors, thresh_sw_str; parse("DENSITY_THRESHOLD",thresh_sw_str);
116 6 : threshold_switch.set(thresh_sw_str, errors );
117 6 : if( errors.length()>0 ) error("errors in DENSITY_THRESHOLD switching function : " + errors );
118 18 : log.printf(" threshold on density of atoms in cylinder equals %s\n",threshold_switch.description().c_str() );
119 :
120 18 : for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
121 18 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
122 12 : double r=cylinder_sw(i,j).get_d0() + cylinder_sw(i,j).get_r0();
123 12 : cell_volume(i,j)=binw_mat(i,j)*pi*r*r;
124 : }
125 : }
126 :
127 : // Find the largest sf cutoff
128 6 : lsfmax=low_sf(0,0).get_dmax();
129 6 : double sfmax=switchingFunction(0,0).get_dmax();
130 6 : double rfmax=cylinder_sw(0,0).get_dmax();
131 18 : for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
132 18 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
133 6 : double tsf=switchingFunction(i,j).get_dmax();
134 6 : if( tsf>sfmax ) sfmax=tsf;
135 6 : double rsf=cylinder_sw(i,j).get_dmax();
136 6 : if( rsf>rfmax ) rfmax=rsf;
137 6 : double lsf=low_sf(i,j).get_dmax();
138 6 : if( lsf>lsfmax ) lsfmax=lsf;
139 : }
140 : }
141 : // Get the width of the bead
142 6 : HistogramBead bead; bead.isNotPeriodic();
143 6 : bead.setKernelType( kerneltype ); bead.set( 0.0, 1.0, sigma );
144 6 : beadrad = bead.getCutoff();
145 :
146 : // Set the link cell cutoff
147 6 : log.printf(" setting cutoff %f \n", sfmax );
148 6 : setLinkCellCutoff( sfmax, std::numeric_limits<double>::max() );
149 :
150 : double maxsize=0;
151 18 : for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
152 18 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
153 6 : if( binw_mat(i,j)>maxsize ) maxsize=binw_mat(i,j);
154 : }
155 : }
156 : // Set the maximum number of bins that we will need to compute
157 6 : maxbins = std::floor( sfmax / maxsize ) + 1;
158 : // Need to resize functions again here to ensure that vector sizes
159 : // are set correctly in AdjacencyMatrixVessel
160 6 : resizeFunctions();
161 6 : }
162 :
163 100 : unsigned TopologyMatrix::getNumberOfQuantities() const {
164 100 : return maxbins+3;
165 : }
166 :
167 24 : void TopologyMatrix::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
168 24 : plumed_assert( id<4 );
169 24 : if( id==0 ) {
170 6 : std::string errors; switchingFunction(j,i).set(desc[0],errors);
171 6 : if( errors.length()!=0 ) error("problem reading switching function description " + errors);
172 6 : if( j!=i) switchingFunction(i,j).set(desc[0],errors);
173 24 : log.printf(" %u th and %u th multicolvar groups must be within %s\n",i+1,j+1,(switchingFunction(i,j).description()).c_str() );
174 18 : } else if( id==1 ) {
175 6 : std::string errors; cylinder_sw(j,i).set(desc[0],errors);
176 6 : if( errors.length()!=0 ) error("problem reading switching function description " + errors);
177 6 : if( j!=i) cylinder_sw(i,j).set(desc[0],errors);
178 24 : log.printf(" there must be not atoms within the cylinder connections atoms of multicolvar groups %u th and %u th. This cylinder has radius %s \n",i+1,j+1,(cylinder_sw(i,j).description()).c_str() );
179 12 : } else if( id==2 ) {
180 6 : std::string errors; low_sf(j,i).set(desc[0],errors);
181 6 : if( errors.length()!=0 ) error("problem reading switching function description " + errors);
182 6 : if( j!=i ) low_sf(i,j).set(desc[0],errors);
183 24 : log.printf(" %u th and %u th multicolvar groups must be further apart than %s\n",i+1,j+1,(low_sf(j,i).description()).c_str() );
184 6 : } else if( id==3 ) {
185 6 : Tools::convert( desc[0], binw_mat(j,i) );
186 6 : if( i!=j ) binw_mat(i,j)=binw_mat(j,i);
187 12 : log.printf(" cylinder for %u th and %u th multicolvar groups is split into bins of length %f \n",i,j,binw_mat(i,j) );
188 : }
189 24 : }
190 :
191 31203 : double TopologyMatrix::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
192 62406 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
193 62406 : if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax2() ) return 1.0;
194 21587 : return 0.0;
195 : }
196 :
197 9616 : double TopologyMatrix::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
198 19232 : HistogramBead bead; bead.isNotPeriodic(); bead.setKernelType( kerneltype );
199 :
200 : // Initialise to zero density on all bins
201 86784 : for(unsigned bin=0; bin<maxbins; ++bin) myatoms.setValue(bin+1,0);
202 : // Calculate whether or not atoms 1 and 2 are within cutoff (can use delta here as pbc are done in atom setup)
203 19232 : Vector d1 = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) ); double d1_len = d1.modulo();
204 9616 : d1 = d1 / d1_len; // Convert vector into director
205 9616 : AtomNumber a1 = myatoms.getAbsoluteIndex( 0 );
206 9616 : AtomNumber a2 = myatoms.getAbsoluteIndex( 1 );
207 235854830 : for(unsigned i=2; i<myatoms.getNumberOfAtoms(); ++i) {
208 117917799 : AtomNumber a3 = myatoms.getAbsoluteIndex( i );
209 235835598 : if( a3!=a1 && a3!=a2 ) calculateForThreeAtoms( i, d1, d1_len, bead, myatoms );
210 : }
211 : // std::vector<double> binvals( 1+maxbins ); for(unsigned i=1;i<maxbins;++i) binvals[i]=myatoms.getValue(i);
212 : // unsigned ii; double fdf;
213 : //std::cout<<"HELLO DENSITY "<<myatoms.getIndex(0)<<" "<<myatoms.getIndex(1)<<" "<<transformStoredValues( binvals, ii, fdf )<<std::endl;
214 :
215 : // Now find the element for which the density is maximal
216 : unsigned vout=2; double max=myatoms.getValue( 2 );
217 67552 : for(unsigned i=3; i<myatoms.getUnderlyingMultiValue().getNumberOfValues()-1; ++i) {
218 28968 : if( myatoms.getValue(i)>max ) { max=myatoms.getValue(i); vout=i; }
219 : }
220 : // Calculate value and derivative of switching function between atoms 1 and 2
221 9616 : double dfuncl, sw = switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ),
222 9616 : getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( d1_len, dfuncl );
223 : // Transform the density
224 9616 : double df, tsw = threshold_switch.calculate( max, df );
225 9616 : if( !doNotCalculateDerivatives() ) {
226 : // Factor of d1_len is required here because d1 is normalized
227 60 : d1 *= d1_len;
228 60 : addAtomDerivatives( 2+maxbins, 0, -dfuncl*d1, myatoms );
229 60 : addAtomDerivatives( 2+maxbins, 1, dfuncl*d1, myatoms );
230 60 : myatoms.addBoxDerivatives( 2+maxbins, (-dfuncl)*Tensor(d1,d1) );
231 : // Update active atoms so that next bit works
232 60 : updateActiveAtoms( myatoms );
233 : // Now finish caclulation of derivatives
234 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
235 25998 : for(unsigned jd=0; jd<myvals.getNumberActive(); ++jd) {
236 12969 : unsigned ider=myvals.getActiveIndex(jd);
237 38907 : myvals.addDerivative( 1, ider, sw*df*max*myvals.getDerivative( vout, ider ) + tsw*myvals.getDerivative( 2+maxbins, ider ) );
238 : }
239 : }
240 9616 : return sw*tsw;
241 : }
242 :
243 117917799 : void TopologyMatrix::calculateForThreeAtoms( const unsigned& iat, const Vector& d1, const double& d1_len,
244 : HistogramBead& bead, multicolvar::AtomValuePack& myatoms ) const {
245 : // Calculate if there are atoms in the cylinder (can use delta here as pbc are done in atom setup)
246 235835598 : Vector d2 = getSeparation( myatoms.getPosition(0), myatoms.getPosition(iat) );
247 : // Now calculate projection of d2 on d1
248 117917799 : double proj=dotProduct(d2,d1);
249 : // This tells us if we are outside the end of the cylinder
250 117917799 : double excess = proj - d1_len;
251 : // Return if we are outside of the cylinder as calculated based on excess
252 212090079 : if( excess>low_sf( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax() ) return;
253 : // Find the length of the cylinder
254 : double binw = binw_mat( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) );
255 81802351 : double lcylinder = (std::floor( d1_len / binw ) + 1)*binw;
256 : // Return if the projection is outside the length of interest
257 81802351 : if( proj<-bead.getCutoff() || proj>(lcylinder+bead.getCutoff()) ) return;
258 :
259 : // Calculate the excess swiching function
260 23745519 : double edf, eval = low_sf( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( excess, edf );
261 : // Calculate the projection on the perpendicular distance from the center of the tube
262 23745519 : double cm = d2.modulo2() - proj*proj;
263 :
264 : // Now calculate the density in the cylinder
265 23745519 : if( cm<cylinder_sw( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax2() ) {
266 104467 : double dfuncr, val = cylinder_sw( getBaseColvarNumber( myatoms.getIndex(0) ),
267 104467 : getBaseColvarNumber( myatoms.getIndex(1) ) ).calculateSqr( cm, dfuncr );
268 : double cellv = cell_volume( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) );
269 104467 : Vector dc1, dc2, dc3, dd1, dd2, dd3, de1, de2, de3;
270 104467 : if( !doNotCalculateDerivatives() ) {
271 4023 : Tensor d1_a1;
272 : // Derivative of director connecting atom1 - atom2 wrt the position of atom 1
273 4023 : d1_a1(0,0) = ( -(d1[1]*d1[1]+d1[2]*d1[2])/d1_len ); // dx/dx
274 4023 : d1_a1(0,1) = ( d1[0]*d1[1]/d1_len ); // dx/dy
275 4023 : d1_a1(0,2) = ( d1[0]*d1[2]/d1_len ); // dx/dz
276 4023 : d1_a1(1,0) = ( d1[1]*d1[0]/d1_len ); // dy/dx
277 4023 : d1_a1(1,1) = ( -(d1[0]*d1[0]+d1[2]*d1[2])/d1_len ); // dy/dy
278 4023 : d1_a1(1,2) = ( d1[1]*d1[2]/d1_len );
279 4023 : d1_a1(2,0) = ( d1[2]*d1[0]/d1_len );
280 4023 : d1_a1(2,1) = ( d1[2]*d1[1]/d1_len );
281 4023 : d1_a1(2,2) = ( -(d1[1]*d1[1]+d1[0]*d1[0])/d1_len );
282 :
283 : // Calculate derivatives of dot product
284 4023 : dd1 = matmul(d2, d1_a1) - d1;
285 4023 : dd2 = matmul(d2, -d1_a1);
286 4023 : dd3 = d1;
287 :
288 : // Calculate derivatives of cross product
289 4023 : dc1 = dfuncr*( -d2 - proj*dd1 );
290 4023 : dc2 = dfuncr*( -proj*dd2 );
291 4023 : dc3 = dfuncr*( d2 - proj*dd3 );
292 :
293 : // Calculate derivatives of excess
294 4023 : de1 = edf*excess*( dd1 + d1 );
295 4023 : de2 = edf*excess*( dd2 - d1 );
296 4023 : de3 = edf*excess*dd3;
297 : }
298 :
299 208934 : Vector pos1 = myatoms.getPosition(0) + d1_len*d1;
300 104467 : Vector pos2 = myatoms.getPosition(0) + d2;
301 104467 : Vector g1derivf,g2derivf,lderivf; Tensor vir;
302 530381 : for(unsigned bin=0; bin<maxbins; ++bin) {
303 425914 : bead.set( bin*binw, (bin+1)*binw, sigma );
304 719410 : if( proj<(bin*binw-bead.getCutoff()) || proj>binw*(bin+1)+bead.getCutoff() ) continue;
305 132418 : double der, contr=bead.calculateWithCutoff( proj, der ) / cellv; der /= cellv;
306 132418 : myatoms.addValue( 2+bin, contr*val*eval );
307 :
308 132418 : if( !doNotCalculateDerivatives() ) {
309 5835 : g1derivf=contr*eval*dc1 + val*eval*der*dd1 + contr*val*de1;
310 5835 : addAtomDerivatives( 2+bin, 0, g1derivf, myatoms );
311 5835 : g2derivf=contr*eval*dc2 + val*eval*der*dd2 + contr*val*de2;
312 5835 : addAtomDerivatives( 2+bin, 1, g2derivf, myatoms );
313 5835 : lderivf=contr*eval*dc3 + val*eval*der*dd3 + contr*val*de3;
314 5835 : addAtomDerivatives( 2+bin, iat, lderivf, myatoms );
315 : // Virial
316 11670 : vir = -Tensor( myatoms.getPosition(0), g1derivf ) - Tensor( pos1, g2derivf ) - Tensor( pos2, lderivf );
317 5835 : myatoms.addBoxDerivatives( 2+bin, vir );
318 : }
319 : }
320 : }
321 : }
322 :
323 : }
324 4839 : }
325 :
|
