Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "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 : 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 override;
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 ) override;
69 : /// This actually calculates the value of the contact function
70 : double calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const override;
71 : /// This does nothing
72 : double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const override;
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 10431 : PLUMED_REGISTER_ACTION(TopologyMatrix,"TOPOLOGY_MATRIX")
79 :
80 7 : void TopologyMatrix::registerKeywords( Keywords& keys ) {
81 7 : AdjacencyMatrixBase::registerKeywords( keys );
82 14 : 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 14 : keys.add("atoms","ATOMS","");
88 14 : keys.add("numbered","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
89 : "The following provides information on the \\ref switchingfunction that are available.");
90 14 : keys.add("numbered","RADIUS","");
91 14 : keys.add("numbered","CYLINDER_SWITCH","a switching function on \\f$(r_{ij}\\cdot r_{ik}-1)/r_{ij}\\f$");
92 14 : keys.add("numbered","BIN_SIZE","the size to use for the bins");
93 21 : keys.add("compulsory","DENSITY_THRESHOLD","");
94 14 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
95 14 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
96 14 : keys.add("hidden","FAKE","");
97 7 : }
98 :
99 6 : TopologyMatrix::TopologyMatrix( const ActionOptions& ao ):
100 : Action(ao),
101 : AdjacencyMatrixBase(ao),
102 12 : maxbins(0)
103 : {
104 12 : readMaxThreeSpeciesMatrix("NODES", "FAKE", "FAKE", "ATOMS", true );
105 6 : unsigned nrows, ncols, ndonor_types; retrieveTypeDimensions( nrows, ncols, ndonor_types );
106 6 : switchingFunction.resize( nrows, ncols ); parseConnectionDescriptions("SWITCH",false,ndonor_types);
107 6 : cylinder_sw.resize( nrows, ncols ); parseConnectionDescriptions("RADIUS",false,ndonor_types);
108 6 : low_sf.resize( nrows, ncols ); parseConnectionDescriptions("CYLINDER_SWITCH",false,ndonor_types);
109 6 : binw_mat.resize( nrows, ncols ); cell_volume.resize( nrows, ncols );
110 6 : parseConnectionDescriptions("BIN_SIZE",false,ndonor_types);
111 : // Read in stuff for grid
112 18 : parse("SIGMA",sigma); parse("KERNEL",kerneltype);
113 : // Read in threshold for density cutoff
114 6 : std::string errors, thresh_sw_str; parse("DENSITY_THRESHOLD",thresh_sw_str);
115 6 : threshold_switch.set(thresh_sw_str, errors );
116 6 : if( errors.length()>0 ) error("errors in DENSITY_THRESHOLD switching function : " + errors );
117 6 : log.printf(" threshold on density of atoms in cylinder equals %s\n",threshold_switch.description().c_str() );
118 :
119 12 : for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
120 12 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
121 6 : double r=cylinder_sw(i,j).get_d0() + cylinder_sw(i,j).get_r0();
122 6 : cell_volume(i,j)=binw_mat(i,j)*pi*r*r;
123 : }
124 : }
125 :
126 : // Find the largest sf cutoff
127 6 : lsfmax=low_sf(0,0).get_dmax();
128 6 : double sfmax=switchingFunction(0,0).get_dmax();
129 6 : double rfmax=cylinder_sw(0,0).get_dmax();
130 12 : for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
131 12 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
132 6 : double tsf=switchingFunction(i,j).get_dmax();
133 6 : if( tsf>sfmax ) sfmax=tsf;
134 6 : double rsf=cylinder_sw(i,j).get_dmax();
135 6 : if( rsf>rfmax ) rfmax=rsf;
136 6 : double lsf=low_sf(i,j).get_dmax();
137 6 : if( lsf>lsfmax ) lsfmax=lsf;
138 : }
139 : }
140 : // Get the width of the bead
141 6 : HistogramBead bead; bead.isNotPeriodic();
142 6 : bead.setKernelType( kerneltype ); bead.set( 0.0, 1.0, sigma );
143 6 : beadrad = bead.getCutoff();
144 :
145 : // Set the link cell cutoff
146 6 : log.printf(" setting cutoff %f \n", sfmax );
147 6 : setLinkCellCutoff( sfmax, std::numeric_limits<double>::max() );
148 :
149 : double maxsize=0;
150 12 : for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
151 12 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
152 6 : if( binw_mat(i,j)>maxsize ) maxsize=binw_mat(i,j);
153 : }
154 : }
155 : // Set the maximum number of bins that we will need to compute
156 6 : maxbins = std::floor( sfmax / maxsize ) + 1;
157 : // Need to resize functions again here to ensure that vector sizes
158 : // are set correctly in AdjacencyMatrixVessel
159 6 : resizeFunctions();
160 6 : }
161 :
162 100 : unsigned TopologyMatrix::getNumberOfQuantities() const {
163 100 : return maxbins+3;
164 : }
165 :
166 24 : void TopologyMatrix::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
167 24 : plumed_assert( id<4 );
168 24 : if( id==0 ) {
169 6 : std::string errors; switchingFunction(j,i).set(desc[0],errors);
170 6 : if( errors.length()!=0 ) error("problem reading switching function description " + errors);
171 6 : if( j!=i) switchingFunction(i,j).set(desc[0],errors);
172 12 : log.printf(" %u th and %u th multicolvar groups must be within %s\n",i+1,j+1,(switchingFunction(i,j).description()).c_str() );
173 18 : } else if( id==1 ) {
174 6 : std::string errors; cylinder_sw(j,i).set(desc[0],errors);
175 6 : if( errors.length()!=0 ) error("problem reading switching function description " + errors);
176 6 : if( j!=i) cylinder_sw(i,j).set(desc[0],errors);
177 12 : 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() );
178 12 : } else if( id==2 ) {
179 6 : std::string errors; low_sf(j,i).set(desc[0],errors);
180 6 : if( errors.length()!=0 ) error("problem reading switching function description " + errors);
181 6 : if( j!=i ) low_sf(i,j).set(desc[0],errors);
182 12 : 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() );
183 6 : } else if( id==3 ) {
184 6 : Tools::convert( desc[0], binw_mat(j,i) );
185 6 : if( i!=j ) binw_mat(i,j)=binw_mat(j,i);
186 6 : 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) );
187 : }
188 24 : }
189 :
190 31203 : double TopologyMatrix::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
191 31203 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
192 62406 : if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax2() ) return 1.0;
193 : return 0.0;
194 : }
195 :
196 9616 : double TopologyMatrix::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
197 9616 : HistogramBead bead; bead.isNotPeriodic(); bead.setKernelType( kerneltype );
198 :
199 : // Initialise to zero density on all bins
200 48200 : for(unsigned bin=0; bin<maxbins; ++bin) myatoms.setValue(bin+1,0);
201 : // Calculate whether or not atoms 1 and 2 are within cutoff (can use delta here as pbc are done in atom setup)
202 9616 : Vector d1 = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) ); double d1_len = d1.modulo();
203 9616 : d1 = d1 / d1_len; // Convert vector into director
204 9616 : AtomNumber a1 = myatoms.getAbsoluteIndex( 0 );
205 9616 : AtomNumber a2 = myatoms.getAbsoluteIndex( 1 );
206 117927415 : for(unsigned i=2; i<myatoms.getNumberOfAtoms(); ++i) {
207 117917799 : AtomNumber a3 = myatoms.getAbsoluteIndex( i );
208 117917799 : if( a3!=a1 && a3!=a2 ) calculateForThreeAtoms( i, d1, d1_len, bead, myatoms );
209 : }
210 : // std::vector<double> binvals( 1+maxbins ); for(unsigned i=1;i<maxbins;++i) binvals[i]=myatoms.getValue(i);
211 : // unsigned ii; double fdf;
212 : //std::cout<<"HELLO DENSITY "<<myatoms.getIndex(0)<<" "<<myatoms.getIndex(1)<<" "<<transformStoredValues( binvals, ii, fdf )<<std::endl;
213 :
214 : // Now find the element for which the density is maximal
215 : unsigned vout=2; double max=myatoms.getValue( 2 );
216 38584 : for(unsigned i=3; i<myatoms.getUnderlyingMultiValue().getNumberOfValues()-1; ++i) {
217 28968 : if( myatoms.getValue(i)>max ) { max=myatoms.getValue(i); vout=i; }
218 : }
219 : // Calculate value and derivative of switching function between atoms 1 and 2
220 : double dfuncl, sw = switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ),
221 9616 : getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( d1_len, dfuncl );
222 : // Transform the density
223 9616 : double df, tsw = threshold_switch.calculate( max, df );
224 9616 : if( !doNotCalculateDerivatives() ) {
225 : // Factor of d1_len is required here because d1 is normalized
226 60 : d1 *= d1_len;
227 60 : addAtomDerivatives( 2+maxbins, 0, -dfuncl*d1, myatoms );
228 60 : addAtomDerivatives( 2+maxbins, 1, dfuncl*d1, myatoms );
229 60 : myatoms.addBoxDerivatives( 2+maxbins, (-dfuncl)*Tensor(d1,d1) );
230 : // Update active atoms so that next bit works
231 60 : updateActiveAtoms( myatoms );
232 : // Now finish caclulation of derivatives
233 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
234 13029 : for(unsigned jd=0; jd<myvals.getNumberActive(); ++jd) {
235 12969 : unsigned ider=myvals.getActiveIndex(jd);
236 12969 : myvals.addDerivative( 1, ider, sw*df*max*myvals.getDerivative( vout, ider ) + tsw*myvals.getDerivative( 2+maxbins, ider ) );
237 : }
238 : }
239 9616 : return sw*tsw;
240 : }
241 :
242 117917799 : void TopologyMatrix::calculateForThreeAtoms( const unsigned& iat, const Vector& d1, const double& d1_len,
243 : HistogramBead& bead, multicolvar::AtomValuePack& myatoms ) const {
244 : // Calculate if there are atoms in the cylinder (can use delta here as pbc are done in atom setup)
245 117917799 : Vector d2 = getSeparation( myatoms.getPosition(0), myatoms.getPosition(iat) );
246 : // Now calculate projection of d2 on d1
247 117917799 : double proj=dotProduct(d2,d1);
248 : // This tells us if we are outside the end of the cylinder
249 117917799 : double excess = proj - d1_len;
250 : // Return if we are outside of the cylinder as calculated based on excess
251 175974631 : if( excess>low_sf( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax() ) return;
252 : // Find the length of the cylinder
253 81802351 : double binw = binw_mat( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) );
254 81802351 : double lcylinder = (std::floor( d1_len / binw ) + 1)*binw;
255 : // Return if the projection is outside the length of interest
256 81802351 : if( proj<-bead.getCutoff() || proj>(lcylinder+bead.getCutoff()) ) return;
257 :
258 : // Calculate the excess swiching function
259 23745519 : double edf, eval = low_sf( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( excess, edf );
260 : // Calculate the projection on the perpendicular distance from the center of the tube
261 23745519 : double cm = d2.modulo2() - proj*proj;
262 :
263 : // Now calculate the density in the cylinder
264 23745519 : if( cm<cylinder_sw( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax2() ) {
265 : double dfuncr, val = cylinder_sw( getBaseColvarNumber( myatoms.getIndex(0) ),
266 104467 : getBaseColvarNumber( myatoms.getIndex(1) ) ).calculateSqr( cm, dfuncr );
267 104467 : double cellv = cell_volume( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) );
268 104467 : Vector dc1, dc2, dc3, dd1, dd2, dd3, de1, de2, de3;
269 104467 : if( !doNotCalculateDerivatives() ) {
270 4023 : Tensor d1_a1;
271 : // Derivative of director connecting atom1 - atom2 wrt the position of atom 1
272 4023 : d1_a1(0,0) = ( -(d1[1]*d1[1]+d1[2]*d1[2])/d1_len ); // dx/dx
273 4023 : d1_a1(0,1) = ( d1[0]*d1[1]/d1_len ); // dx/dy
274 4023 : d1_a1(0,2) = ( d1[0]*d1[2]/d1_len ); // dx/dz
275 4023 : d1_a1(1,0) = ( d1[1]*d1[0]/d1_len ); // dy/dx
276 4023 : d1_a1(1,1) = ( -(d1[0]*d1[0]+d1[2]*d1[2])/d1_len ); // dy/dy
277 4023 : d1_a1(1,2) = ( d1[1]*d1[2]/d1_len );
278 4023 : d1_a1(2,0) = ( d1[2]*d1[0]/d1_len );
279 4023 : d1_a1(2,1) = ( d1[2]*d1[1]/d1_len );
280 4023 : d1_a1(2,2) = ( -(d1[1]*d1[1]+d1[0]*d1[0])/d1_len );
281 :
282 : // Calculate derivatives of dot product
283 4023 : dd1 = matmul(d2, d1_a1) - d1;
284 4023 : dd2 = matmul(d2, -d1_a1);
285 4023 : dd3 = d1;
286 :
287 : // Calculate derivatives of cross product
288 4023 : dc1 = dfuncr*( -d2 - proj*dd1 );
289 4023 : dc2 = dfuncr*( -proj*dd2 );
290 4023 : dc3 = dfuncr*( d2 - proj*dd3 );
291 :
292 : // Calculate derivatives of excess
293 4023 : de1 = edf*excess*( dd1 + d1 );
294 4023 : de2 = edf*excess*( dd2 - d1 );
295 4023 : de3 = edf*excess*dd3;
296 : }
297 :
298 104467 : Vector pos1 = myatoms.getPosition(0) + d1_len*d1;
299 104467 : Vector pos2 = myatoms.getPosition(0) + d2;
300 104467 : Vector g1derivf,g2derivf,lderivf; Tensor vir;
301 530381 : for(unsigned bin=0; bin<maxbins; ++bin) {
302 425914 : bead.set( bin*binw, (bin+1)*binw, sigma );
303 425914 : if( proj<(bin*binw-bead.getCutoff()) || proj>binw*(bin+1)+bead.getCutoff() ) continue;
304 132418 : double der, contr=bead.calculateWithCutoff( proj, der ) / cellv; der /= cellv;
305 132418 : myatoms.addValue( 2+bin, contr*val*eval );
306 :
307 132418 : if( !doNotCalculateDerivatives() ) {
308 5835 : g1derivf=contr*eval*dc1 + val*eval*der*dd1 + contr*val*de1;
309 5835 : addAtomDerivatives( 2+bin, 0, g1derivf, myatoms );
310 5835 : g2derivf=contr*eval*dc2 + val*eval*der*dd2 + contr*val*de2;
311 5835 : addAtomDerivatives( 2+bin, 1, g2derivf, myatoms );
312 5835 : lderivf=contr*eval*dc3 + val*eval*der*dd3 + contr*val*de3;
313 5835 : addAtomDerivatives( 2+bin, iat, lderivf, myatoms );
314 : // Virial
315 5835 : vir = -Tensor( myatoms.getPosition(0), g1derivf ) - Tensor( pos1, g2derivf ) - Tensor( pos2, lderivf );
316 5835 : myatoms.addBoxDerivatives( 2+bin, vir );
317 : }
318 : }
319 : }
320 : }
321 :
322 : }
323 : }
324 :
|
