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 : {
87 16 : std::string sfinput,errors; parse("SWITCH",sfinput);
88 8 : if( sfinput.length()==0 ) error("could not find SWITCH keyword");
89 8 : switchingFunction.set(sfinput,errors);
90 8 : if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
91 :
92 16 : std::string hsfinput; parse("CYLINDER_SWITCH",hsfinput);
93 8 : if( hsfinput.length()==0 ) error("could not find CYLINDER_SWITCH keyword");
94 8 : low_sf.set(hsfinput,errors);
95 8 : if( errors.length()!=0 ) error("problem reading CYLINDER_SWITCH keyword : " + errors );
96 :
97 16 : std::string asfinput; parse("RADIUS",asfinput);
98 8 : if( asfinput.length()==0 ) error("could not find RADIUS keyword");
99 8 : cylinder_sw.set(asfinput,errors);
100 8 : if( errors.length()!=0 ) error("problem reading RADIUS keyword : " + errors );
101 :
102 16 : std::string tsfinput; parse("DENSITY_THRESHOLD",tsfinput);
103 8 : if( tsfinput.length()==0 ) error("could not find DENSITY_THRESHOLD keyword");
104 8 : threshold_switch.set(tsfinput,errors);
105 8 : if( errors.length()!=0 ) error("problem reading DENSITY_THRESHOLD keyword : " + errors );
106 : // Read in stuff for grid
107 24 : parse("SIGMA",sigma); parse("KERNEL",kerneltype); parse("BIN_SIZE",binw_mat);
108 :
109 : // Set the link cell cutoff
110 8 : setLinkCellCutoff( true, switchingFunction.get_dmax(), std::numeric_limits<double>::max() );
111 : // Set the number of bins
112 8 : maxbins = std::floor( switchingFunction.get_dmax() / binw_mat ) + 1;
113 : // Set the cell volume
114 8 : double r=cylinder_sw.get_d0() + cylinder_sw.get_r0();
115 8 : cell_volume=binw_mat*pi*r*r;
116 :
117 : // And check everything has been read in correctly
118 8 : checkRead();
119 8 : }
120 :
121 69150 : double TopologyMatrix::calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const {
122 : // Compute switching function on distance between atoms
123 69150 : Vector distance = pbcDistance( pos1, pos2 ); double len2 = distance.modulo2();
124 69150 : if( len2>switchingFunction.get_dmax2() ) return 0.0;
125 26332 : double dfuncl, sw = switchingFunction.calculateSqr( len2, dfuncl );
126 :
127 : // Now run through all sea atoms
128 26332 : HistogramBead bead; bead.isNotPeriodic(); bead.setKernelType( kerneltype );
129 26332 : Vector g1derivf,g2derivf,lderivf; Tensor vir; double binlength = maxbins * binw_mat;
130 26332 : MultiValue tvals( maxbins, myvals.getNumberOfDerivatives() );
131 236742358 : for(unsigned i=0; i<natoms; ++i) {
132 : // Position of sea atom (this will be the origin)
133 : Vector d2 = getPosition(i,myvals);
134 : // Vector connecting sea atom and first in bond taking pbc into account
135 : Vector d20 = pbcDistance( d2, pos1 );
136 : // Vector connecting sea atom and second in bond taking pbc into account
137 : Vector d21 = pbcDistance( d2, pos2 );
138 : // Now length of bond modulus and so on -- no pbc here as we want sea atom in middle
139 236716026 : Vector d1 = delta( d20, d21 ); double d1_len = d1.modulo(); d1 = d1 / d1_len;
140 : // Switching function on distance between nodes
141 371343234 : if( d1_len>switchingFunction.get_dmax() ) continue ;
142 : // Ensure that the center of the bins are on the center of the bond connecting the two atoms
143 183588938 : double start2atom = 0.5*(binlength-d1_len); Vector dstart = d20 - start2atom*d1;
144 : // Now calculate projection of axis of cylinder
145 183588938 : double proj=dotProduct(-dstart,d1);
146 : // Calculate length of vector connecting start of cylinder to first atom
147 : // And now work out projection on vector connecting start and end of cylinder
148 183588938 : double proj_between = proj - start2atom;
149 : // This tells us if we are outside the end of the cylinder
150 183588938 : double excess = proj_between - d1_len;
151 : // Return if we are outside of the cylinder as calculated based on excess
152 183588938 : if( excess>low_sf.get_dmax() || -proj_between>low_sf.get_dmax() ) continue;
153 : // Calculate the excess swiching functions
154 48961730 : double edf1, eval1 = low_sf.calculate( excess, edf1 );
155 48961730 : double edf2, eval2 = low_sf.calculate( -proj_between, edf2 );
156 : // Calculate the projection on the perpendicular distance from the center of the tube
157 48961730 : double cm = dstart.modulo2() - proj*proj;
158 :
159 : // Now calculate the density in the cylinder
160 48961730 : if( cm>0 && cm<cylinder_sw.get_dmax2() ) {
161 288350 : double dfuncr, val = cylinder_sw.calculateSqr( cm, dfuncr );
162 288350 : Vector dc1, dc2, dc3, dd1, dd2, dd3, de1, de2, de3;
163 288350 : if( !doNotCalculateDerivatives() ) {
164 28284 : Tensor d1_a1;
165 : // Derivative of director connecting atom1 - atom2 wrt the position of atom 1
166 28284 : d1_a1(0,0) = ( -(d1[1]*d1[1]+d1[2]*d1[2])/d1_len ); // dx/dx
167 28284 : d1_a1(0,1) = ( d1[0]*d1[1]/d1_len ); // dx/dy
168 28284 : d1_a1(0,2) = ( d1[0]*d1[2]/d1_len ); // dx/dz
169 28284 : d1_a1(1,0) = ( d1[1]*d1[0]/d1_len ); // dy/dx
170 28284 : d1_a1(1,1) = ( -(d1[0]*d1[0]+d1[2]*d1[2])/d1_len ); // dy/dy
171 28284 : d1_a1(1,2) = ( d1[1]*d1[2]/d1_len );
172 28284 : d1_a1(2,0) = ( d1[2]*d1[0]/d1_len );
173 28284 : d1_a1(2,1) = ( d1[2]*d1[1]/d1_len );
174 28284 : d1_a1(2,2) = ( -(d1[1]*d1[1]+d1[0]*d1[0])/d1_len );
175 :
176 : // Calculate derivatives of dot product
177 28284 : dd1 = matmul(-dstart, d1_a1) - 0.5*d1;
178 28284 : dd2 = matmul(-dstart, -d1_a1) - 0.5*d1;
179 28284 : dd3 = d1;
180 :
181 : // Calculate derivatives of cross product
182 28284 : Vector der( -0.5*binlength*matmul( d1_a1,dstart ) );
183 28284 : dc1 = dfuncr*( 0.5*dstart + der - proj*dd1 );
184 28284 : dc2 = dfuncr*( 0.5*dstart - der - proj*dd2 );
185 28284 : dc3 = dfuncr*( -dstart - proj*dd3 );
186 :
187 : // Calculate derivatives of excess
188 28284 : de1 = eval2*edf1*excess*(dd1 + 0.5*d1 ) + eval1*edf2*proj_between*(dd1 - 0.5*d1);
189 28284 : de2 = eval2*edf1*excess*(dd2 - 0.5*d1 ) + eval1*edf2*proj_between*(dd2 + 0.5*d1);
190 28284 : de3 = ( eval2*edf1*excess + eval1*edf2*proj_between )*dd3;
191 : }
192 2225914 : for(unsigned bin=0; bin<maxbins; ++bin) {
193 1937564 : bead.set( bin*binw_mat, (bin+1)*binw_mat, sigma );
194 1937564 : if( proj<(bin*binw_mat-bead.getCutoff()) || proj>binw_mat*(bin+1)+bead.getCutoff() ) continue;
195 386996 : double der, contr=bead.calculateWithCutoff( proj, der ) / cell_volume; der /= cell_volume;
196 386996 : tvals.addValue( bin, contr*val*eval1*eval2 );
197 :
198 386996 : if( !doNotCalculateDerivatives() ) {
199 38132 : g1derivf=contr*eval1*eval2*dc1 + val*eval1*eval2*der*dd1 + contr*val*de1;
200 38132 : tvals.addDerivative( bin, 3*myvals.getTaskIndex()+0, g1derivf[0] );
201 38132 : tvals.addDerivative( bin, 3*myvals.getTaskIndex()+1, g1derivf[1] );
202 38132 : tvals.addDerivative( bin, 3*myvals.getTaskIndex()+2, g1derivf[2] );
203 38132 : g2derivf=contr*eval1*eval2*dc2 + val*eval1*eval2*der*dd2 + contr*val*de2;
204 38132 : tvals.addDerivative( bin, 3*myvals.getSecondTaskIndex()+0, g2derivf[0] );
205 38132 : tvals.addDerivative( bin, 3*myvals.getSecondTaskIndex()+1, g2derivf[1] );
206 38132 : tvals.addDerivative( bin, 3*myvals.getSecondTaskIndex()+2, g2derivf[2] );
207 38132 : lderivf=contr*eval1*eval2*dc3 + val*eval1*eval2*der*dd3 + contr*val*de3;
208 38132 : unsigned tindex = myvals.getIndices()[ i + myvals.getSplitIndex() ];
209 38132 : tvals.addDerivative( bin, 3*tindex+0, lderivf[0] );
210 38132 : tvals.addDerivative( bin, 3*tindex+1, lderivf[1] );
211 38132 : tvals.addDerivative( bin, 3*tindex+2, lderivf[2] );
212 : // Virial
213 38132 : vir = - Tensor( d20, g1derivf ) - Tensor( d21, g2derivf );
214 38132 : unsigned nbase = 3*getNumberOfAtoms();
215 38132 : tvals.addDerivative( bin, nbase+0, vir(0,0) );
216 38132 : tvals.addDerivative( bin, nbase+1, vir(0,1) );
217 38132 : tvals.addDerivative( bin, nbase+2, vir(0,2) );
218 38132 : tvals.addDerivative( bin, nbase+3, vir(1,0) );
219 38132 : tvals.addDerivative( bin, nbase+4, vir(1,1) );
220 38132 : tvals.addDerivative( bin, nbase+5, vir(1,2) );
221 38132 : tvals.addDerivative( bin, nbase+6, vir(2,0) );
222 38132 : tvals.addDerivative( bin, nbase+7, vir(2,1) );
223 38132 : tvals.addDerivative( bin, nbase+8, vir(2,2) );
224 : }
225 : }
226 : }
227 : }
228 : // Find maximum density
229 : double max = tvals.get(0); unsigned vout = 0;
230 305488 : for(unsigned i=1; i<maxbins; ++i) {
231 279156 : if( tvals.get(i)>max ) { max=tvals.get(i); vout=i; }
232 : }
233 : // Transform the density
234 26332 : double df, tsw = threshold_switch.calculate( max, df );
235 26332 : if( fabs(sw*tsw)<epsilon ) return 0;
236 :
237 3516 : if( !doNotCalculateDerivatives() ) {
238 942 : Vector ader; Tensor vir; Vector ddd = tsw*dfuncl*distance;
239 942 : ader[0] = tvals.getDerivative( vout, 3*myvals.getTaskIndex()+0 );
240 942 : ader[1] = tvals.getDerivative( vout, 3*myvals.getTaskIndex()+1 );
241 942 : ader[2] = tvals.getDerivative( vout, 3*myvals.getTaskIndex()+2 );
242 942 : addAtomDerivatives( 0, sw*df*max*ader - ddd, myvals );
243 942 : ader[0] = tvals.getDerivative( vout, 3*myvals.getSecondTaskIndex()+0 );
244 942 : ader[1] = tvals.getDerivative( vout, 3*myvals.getSecondTaskIndex()+1 );
245 942 : ader[2] = tvals.getDerivative( vout, 3*myvals.getSecondTaskIndex()+2 );
246 942 : addAtomDerivatives( 1, sw*df*max*ader + ddd, myvals );
247 109488 : for(unsigned i=0; i<natoms; ++i) {
248 108546 : unsigned tindex = myvals.getIndices()[ i + myvals.getSplitIndex() ];
249 108546 : ader[0] = tvals.getDerivative( vout, 3*tindex+0 );
250 108546 : ader[1] = tvals.getDerivative( vout, 3*tindex+1 );
251 108546 : ader[2] = tvals.getDerivative( vout, 3*tindex+2 );
252 108546 : addThirdAtomDerivatives( i, sw*df*max*ader, myvals );
253 : }
254 942 : unsigned nbase = 3*getNumberOfAtoms(); Tensor vird(ddd,distance);
255 942 : vir(0,0) = sw*df*max*tvals.getDerivative( vout, nbase+0 ) - vird(0,0);
256 942 : vir(0,1) = sw*df*max*tvals.getDerivative( vout, nbase+1 ) - vird(0,1);
257 942 : vir(0,2) = sw*df*max*tvals.getDerivative( vout, nbase+2 ) - vird(0,2);
258 942 : vir(1,0) = sw*df*max*tvals.getDerivative( vout, nbase+3 ) - vird(1,0);
259 942 : vir(1,1) = sw*df*max*tvals.getDerivative( vout, nbase+4 ) - vird(1,1);
260 942 : vir(1,2) = sw*df*max*tvals.getDerivative( vout, nbase+5 ) - vird(1,2);
261 942 : vir(2,0) = sw*df*max*tvals.getDerivative( vout, nbase+6 ) - vird(2,0);
262 942 : vir(2,1) = sw*df*max*tvals.getDerivative( vout, nbase+7 ) - vird(2,1);
263 942 : vir(2,2) = sw*df*max*tvals.getDerivative( vout, nbase+8 ) - vird(2,2);
264 942 : addBoxDerivatives( vir, myvals );
265 : }
266 : return sw*tsw;
267 26332 : }
268 :
269 : }
270 : }
|
