Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2020 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 "tools/SwitchingFunction.h"
24 : #include "tools/HistogramBead.h"
25 : #include "core/ActionRegister.h"
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : //+PLUMEDOC MATRIX TOPOLOGY_MATRIX
31 : /*
32 : Adjacency matrix in which two atoms are adjacent if they are connected topologically
33 :
34 : \par Examples
35 :
36 :
37 : */
38 : //+ENDPLUMEDOC
39 :
40 : namespace PLMD {
41 : namespace adjmat {
42 :
43 : class TopologyMatrix : public AdjacencyMatrixBase {
44 : private:
45 : /// The width to use for the kernel density estimation and the
46 : /// sizes of the bins to be used in kernel density estimation
47 : double sigma;
48 : std::string kerneltype;
49 : /// The maximum number of bins that will be used
50 : /// This is calculated based on the dmax of the switching functions
51 : unsigned maxbins;
52 : /// The volume of the cells
53 : double cell_volume;
54 : /// switching function
55 : SwitchingFunction switchingFunction;
56 : SwitchingFunction cylinder_sw;
57 : SwitchingFunction low_sf;
58 : double binw_mat;
59 : SwitchingFunction threshold_switch;
60 : public:
61 : static void registerKeywords( Keywords& keys );
62 : explicit TopologyMatrix(const ActionOptions&);
63 : // active methods:
64 : double calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const override;
65 : };
66 :
67 : PLUMED_REGISTER_ACTION(TopologyMatrix,"TOPOLOGY_MATRIX")
68 :
69 10 : void TopologyMatrix::registerKeywords( Keywords& keys ) {
70 10 : AdjacencyMatrixBase::registerKeywords( keys );
71 20 : keys.add("atoms","BACKGROUND_ATOMS","the list of atoms that should be considered as part of the background density");
72 20 : keys.add("compulsory","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
73 : "The following provides information on the \\ref switchingfunction that are available. "
74 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
75 20 : keys.add("compulsory","RADIUS","");
76 20 : keys.add("compulsory","CYLINDER_SWITCH","a switching function on ( r_ij . r_ik - 1 )/r_ij");
77 20 : keys.add("compulsory","BIN_SIZE","the size to use for the bins");
78 20 : keys.add("compulsory","DENSITY_THRESHOLD","");
79 20 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
80 20 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
81 10 : }
82 :
83 8 : TopologyMatrix::TopologyMatrix(const ActionOptions&ao):
84 : Action(ao),
85 8 : AdjacencyMatrixBase(ao) {
86 : std::string sfinput,errors;
87 16 : parse("SWITCH",sfinput);
88 8 : if( sfinput.length()==0 ) {
89 0 : error("could not find SWITCH keyword");
90 : }
91 8 : switchingFunction.set(sfinput,errors);
92 8 : if( errors.length()!=0 ) {
93 0 : error("problem reading SWITCH keyword : " + errors );
94 : }
95 :
96 : std::string hsfinput;
97 16 : parse("CYLINDER_SWITCH",hsfinput);
98 8 : if( hsfinput.length()==0 ) {
99 0 : error("could not find CYLINDER_SWITCH keyword");
100 : }
101 8 : low_sf.set(hsfinput,errors);
102 8 : if( errors.length()!=0 ) {
103 0 : error("problem reading CYLINDER_SWITCH keyword : " + errors );
104 : }
105 :
106 : std::string asfinput;
107 16 : parse("RADIUS",asfinput);
108 8 : if( asfinput.length()==0 ) {
109 0 : error("could not find RADIUS keyword");
110 : }
111 8 : cylinder_sw.set(asfinput,errors);
112 8 : if( errors.length()!=0 ) {
113 0 : error("problem reading RADIUS keyword : " + errors );
114 : }
115 :
116 : std::string tsfinput;
117 16 : parse("DENSITY_THRESHOLD",tsfinput);
118 8 : if( tsfinput.length()==0 ) {
119 0 : error("could not find DENSITY_THRESHOLD keyword");
120 : }
121 8 : threshold_switch.set(tsfinput,errors);
122 8 : if( errors.length()!=0 ) {
123 0 : error("problem reading DENSITY_THRESHOLD keyword : " + errors );
124 : }
125 : // Read in stuff for grid
126 8 : parse("SIGMA",sigma);
127 8 : parse("KERNEL",kerneltype);
128 8 : parse("BIN_SIZE",binw_mat);
129 :
130 : // Set the link cell cutoff
131 8 : setLinkCellCutoff( true, switchingFunction.get_dmax(), std::numeric_limits<double>::max() );
132 : // Set the number of bins
133 8 : maxbins = std::floor( switchingFunction.get_dmax() / binw_mat ) + 1;
134 : // Set the cell volume
135 8 : double r=cylinder_sw.get_d0() + cylinder_sw.get_r0();
136 8 : cell_volume=binw_mat*pi*r*r;
137 :
138 : // And check everything has been read in correctly
139 8 : checkRead();
140 8 : }
141 :
142 69150 : double TopologyMatrix::calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const {
143 : // Compute switching function on distance between atoms
144 : Vector distance = pbcDistance( pos1, pos2 );
145 69150 : double len2 = distance.modulo2();
146 69150 : if( len2>switchingFunction.get_dmax2() ) {
147 : return 0.0;
148 : }
149 26332 : double dfuncl, sw = switchingFunction.calculateSqr( len2, dfuncl );
150 :
151 : // Now run through all sea atoms
152 26332 : HistogramBead bead;
153 : bead.isNotPeriodic();
154 26332 : bead.setKernelType( kerneltype );
155 26332 : Vector g1derivf,g2derivf,lderivf;
156 26332 : Tensor vir;
157 26332 : double binlength = maxbins * binw_mat;
158 26332 : MultiValue tvals( maxbins, myvals.getNumberOfDerivatives() );
159 236742358 : for(unsigned i=0; i<natoms; ++i) {
160 : // Position of sea atom (this will be the origin)
161 : Vector d2 = getPosition(i,myvals);
162 : // Vector connecting sea atom and first in bond taking pbc into account
163 : Vector d20 = pbcDistance( d2, pos1 );
164 : // Vector connecting sea atom and second in bond taking pbc into account
165 : Vector d21 = pbcDistance( d2, pos2 );
166 : // Now length of bond modulus and so on -- no pbc here as we want sea atom in middle
167 236716026 : Vector d1 = delta( d20, d21 );
168 236716026 : double d1_len = d1.modulo();
169 236716026 : d1 = d1 / d1_len;
170 : // Switching function on distance between nodes
171 236716026 : if( d1_len>switchingFunction.get_dmax() ) {
172 187754296 : continue ;
173 : }
174 : // Ensure that the center of the bins are on the center of the bond connecting the two atoms
175 183588938 : double start2atom = 0.5*(binlength-d1_len);
176 183588938 : Vector dstart = d20 - start2atom*d1;
177 : // Now calculate projection of axis of cylinder
178 183588938 : double proj=dotProduct(-dstart,d1);
179 : // Calculate length of vector connecting start of cylinder to first atom
180 : // And now work out projection on vector connecting start and end of cylinder
181 183588938 : double proj_between = proj - start2atom;
182 : // This tells us if we are outside the end of the cylinder
183 183588938 : double excess = proj_between - d1_len;
184 : // Return if we are outside of the cylinder as calculated based on excess
185 183588938 : if( excess>low_sf.get_dmax() || -proj_between>low_sf.get_dmax() ) {
186 134627208 : continue;
187 : }
188 : // Calculate the excess swiching functions
189 48961730 : double edf1, eval1 = low_sf.calculate( excess, edf1 );
190 48961730 : double edf2, eval2 = low_sf.calculate( -proj_between, edf2 );
191 : // Calculate the projection on the perpendicular distance from the center of the tube
192 48961730 : double cm = dstart.modulo2() - proj*proj;
193 :
194 : // Now calculate the density in the cylinder
195 48961730 : if( cm>0 && cm<cylinder_sw.get_dmax2() ) {
196 288350 : double dfuncr, val = cylinder_sw.calculateSqr( cm, dfuncr );
197 288350 : Vector dc1, dc2, dc3, dd1, dd2, dd3, de1, de2, de3;
198 288350 : if( !doNotCalculateDerivatives() ) {
199 28284 : Tensor d1_a1;
200 : // Derivative of director connecting atom1 - atom2 wrt the position of atom 1
201 28284 : d1_a1(0,0) = ( -(d1[1]*d1[1]+d1[2]*d1[2])/d1_len ); // dx/dx
202 28284 : d1_a1(0,1) = ( d1[0]*d1[1]/d1_len ); // dx/dy
203 28284 : d1_a1(0,2) = ( d1[0]*d1[2]/d1_len ); // dx/dz
204 28284 : d1_a1(1,0) = ( d1[1]*d1[0]/d1_len ); // dy/dx
205 28284 : d1_a1(1,1) = ( -(d1[0]*d1[0]+d1[2]*d1[2])/d1_len ); // dy/dy
206 28284 : d1_a1(1,2) = ( d1[1]*d1[2]/d1_len );
207 28284 : d1_a1(2,0) = ( d1[2]*d1[0]/d1_len );
208 28284 : d1_a1(2,1) = ( d1[2]*d1[1]/d1_len );
209 28284 : d1_a1(2,2) = ( -(d1[1]*d1[1]+d1[0]*d1[0])/d1_len );
210 :
211 : // Calculate derivatives of dot product
212 28284 : dd1 = matmul(-dstart, d1_a1) - 0.5*d1;
213 28284 : dd2 = matmul(-dstart, -d1_a1) - 0.5*d1;
214 28284 : dd3 = d1;
215 :
216 : // Calculate derivatives of cross product
217 28284 : Vector der( -0.5*binlength*matmul( d1_a1,dstart ) );
218 28284 : dc1 = dfuncr*( 0.5*dstart + der - proj*dd1 );
219 28284 : dc2 = dfuncr*( 0.5*dstart - der - proj*dd2 );
220 28284 : dc3 = dfuncr*( -dstart - proj*dd3 );
221 :
222 : // Calculate derivatives of excess
223 28284 : de1 = eval2*edf1*excess*(dd1 + 0.5*d1 ) + eval1*edf2*proj_between*(dd1 - 0.5*d1);
224 28284 : de2 = eval2*edf1*excess*(dd2 - 0.5*d1 ) + eval1*edf2*proj_between*(dd2 + 0.5*d1);
225 28284 : de3 = ( eval2*edf1*excess + eval1*edf2*proj_between )*dd3;
226 : }
227 2225914 : for(unsigned bin=0; bin<maxbins; ++bin) {
228 1937564 : bead.set( bin*binw_mat, (bin+1)*binw_mat, sigma );
229 1937564 : if( proj<(bin*binw_mat-bead.getCutoff()) || proj>binw_mat*(bin+1)+bead.getCutoff() ) {
230 1550568 : continue;
231 : }
232 386996 : double der, contr=bead.calculateWithCutoff( proj, der ) / cell_volume;
233 386996 : der /= cell_volume;
234 386996 : tvals.addValue( bin, contr*val*eval1*eval2 );
235 :
236 386996 : if( !doNotCalculateDerivatives() ) {
237 38132 : g1derivf=contr*eval1*eval2*dc1 + val*eval1*eval2*der*dd1 + contr*val*de1;
238 38132 : tvals.addDerivative( bin, 3*myvals.getTaskIndex()+0, g1derivf[0] );
239 38132 : tvals.addDerivative( bin, 3*myvals.getTaskIndex()+1, g1derivf[1] );
240 38132 : tvals.addDerivative( bin, 3*myvals.getTaskIndex()+2, g1derivf[2] );
241 38132 : g2derivf=contr*eval1*eval2*dc2 + val*eval1*eval2*der*dd2 + contr*val*de2;
242 38132 : tvals.addDerivative( bin, 3*myvals.getSecondTaskIndex()+0, g2derivf[0] );
243 38132 : tvals.addDerivative( bin, 3*myvals.getSecondTaskIndex()+1, g2derivf[1] );
244 38132 : tvals.addDerivative( bin, 3*myvals.getSecondTaskIndex()+2, g2derivf[2] );
245 38132 : lderivf=contr*eval1*eval2*dc3 + val*eval1*eval2*der*dd3 + contr*val*de3;
246 38132 : unsigned tindex = myvals.getIndices()[ i + myvals.getSplitIndex() ];
247 38132 : tvals.addDerivative( bin, 3*tindex+0, lderivf[0] );
248 38132 : tvals.addDerivative( bin, 3*tindex+1, lderivf[1] );
249 38132 : tvals.addDerivative( bin, 3*tindex+2, lderivf[2] );
250 : // Virial
251 38132 : vir = - Tensor( d20, g1derivf ) - Tensor( d21, g2derivf );
252 38132 : unsigned nbase = 3*getNumberOfAtoms();
253 38132 : tvals.addDerivative( bin, nbase+0, vir(0,0) );
254 38132 : tvals.addDerivative( bin, nbase+1, vir(0,1) );
255 38132 : tvals.addDerivative( bin, nbase+2, vir(0,2) );
256 38132 : tvals.addDerivative( bin, nbase+3, vir(1,0) );
257 38132 : tvals.addDerivative( bin, nbase+4, vir(1,1) );
258 38132 : tvals.addDerivative( bin, nbase+5, vir(1,2) );
259 38132 : tvals.addDerivative( bin, nbase+6, vir(2,0) );
260 38132 : tvals.addDerivative( bin, nbase+7, vir(2,1) );
261 38132 : tvals.addDerivative( bin, nbase+8, vir(2,2) );
262 : }
263 : }
264 : }
265 : }
266 : // Find maximum density
267 : double max = tvals.get(0);
268 : unsigned vout = 0;
269 305488 : for(unsigned i=1; i<maxbins; ++i) {
270 279156 : if( tvals.get(i)>max ) {
271 : max=tvals.get(i);
272 : vout=i;
273 : }
274 : }
275 : // Transform the density
276 26332 : double df, tsw = threshold_switch.calculate( max, df );
277 26332 : if( fabs(sw*tsw)<epsilon ) {
278 : return 0;
279 : }
280 :
281 3516 : if( !doNotCalculateDerivatives() ) {
282 942 : Vector ader;
283 942 : Tensor vir;
284 942 : Vector ddd = tsw*dfuncl*distance;
285 942 : ader[0] = tvals.getDerivative( vout, 3*myvals.getTaskIndex()+0 );
286 942 : ader[1] = tvals.getDerivative( vout, 3*myvals.getTaskIndex()+1 );
287 942 : ader[2] = tvals.getDerivative( vout, 3*myvals.getTaskIndex()+2 );
288 942 : addAtomDerivatives( 0, sw*df*max*ader - ddd, myvals );
289 942 : ader[0] = tvals.getDerivative( vout, 3*myvals.getSecondTaskIndex()+0 );
290 942 : ader[1] = tvals.getDerivative( vout, 3*myvals.getSecondTaskIndex()+1 );
291 942 : ader[2] = tvals.getDerivative( vout, 3*myvals.getSecondTaskIndex()+2 );
292 942 : addAtomDerivatives( 1, sw*df*max*ader + ddd, myvals );
293 109488 : for(unsigned i=0; i<natoms; ++i) {
294 108546 : unsigned tindex = myvals.getIndices()[ i + myvals.getSplitIndex() ];
295 108546 : ader[0] = tvals.getDerivative( vout, 3*tindex+0 );
296 108546 : ader[1] = tvals.getDerivative( vout, 3*tindex+1 );
297 108546 : ader[2] = tvals.getDerivative( vout, 3*tindex+2 );
298 108546 : addThirdAtomDerivatives( i, sw*df*max*ader, myvals );
299 : }
300 942 : unsigned nbase = 3*getNumberOfAtoms();
301 942 : Tensor vird(ddd,distance);
302 942 : vir(0,0) = sw*df*max*tvals.getDerivative( vout, nbase+0 ) - vird(0,0);
303 942 : vir(0,1) = sw*df*max*tvals.getDerivative( vout, nbase+1 ) - vird(0,1);
304 942 : vir(0,2) = sw*df*max*tvals.getDerivative( vout, nbase+2 ) - vird(0,2);
305 942 : vir(1,0) = sw*df*max*tvals.getDerivative( vout, nbase+3 ) - vird(1,0);
306 942 : vir(1,1) = sw*df*max*tvals.getDerivative( vout, nbase+4 ) - vird(1,1);
307 942 : vir(1,2) = sw*df*max*tvals.getDerivative( vout, nbase+5 ) - vird(1,2);
308 942 : vir(2,0) = sw*df*max*tvals.getDerivative( vout, nbase+6 ) - vird(2,0);
309 942 : vir(2,1) = sw*df*max*tvals.getDerivative( vout, nbase+7 ) - vird(2,1);
310 942 : vir(2,2) = sw*df*max*tvals.getDerivative( vout, nbase+8 ) - vird(2,2);
311 942 : addBoxDerivatives( vir, myvals );
312 : }
313 : return sw*tsw;
314 26332 : }
315 :
316 : }
317 : }
|
