Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "core/ActionRegister.h"
25 :
26 : #include <string>
27 : #include <cmath>
28 :
29 : namespace PLMD {
30 : namespace adjmat {
31 :
32 : //+PLUMEDOC MCOLVAR BRIDGE_MATRIX
33 : /*
34 : Calculate the number of atoms that bridge two parts of a structure
35 :
36 : This quantity calculates:
37 :
38 : \f[
39 : f(x) = \sum_{ijk} s_A(r_{ij})s_B(r_{ik})
40 : \f]
41 :
42 : where the sum over \f$i\f$ is over all the ``bridging atoms" and
43 : \f$s_A\f$ and \f$s_B\f$ are \ref switchingfunction.
44 :
45 : \par Examples
46 :
47 : The following example instructs plumed to calculate the number of water molecules
48 : that are bridging between atoms 1-10 and atoms 11-20 and to print the value
49 : to a file
50 :
51 : \plumedfile
52 : w1: BRIDGE BRIDGING_ATOMS=100-200 GROUPA=1-10 GROUPB=11-20 SWITCH={RATIONAL R_0=0.2}
53 : PRINT ARG=w1 FILE=colvar
54 : \endplumedfile
55 :
56 : */
57 : //+ENDPLUMEDOC
58 :
59 : class BridgeMatrix : public AdjacencyMatrixBase {
60 : private:
61 : Vector dij, dik;
62 : SwitchingFunction sf1;
63 : SwitchingFunction sf2;
64 : public:
65 : static void registerKeywords( Keywords& keys );
66 : explicit BridgeMatrix(const ActionOptions&);
67 : // active methods:
68 : double calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const override;
69 : };
70 :
71 : PLUMED_REGISTER_ACTION(BridgeMatrix,"BRIDGE_MATRIX")
72 :
73 18 : void BridgeMatrix::registerKeywords( Keywords& keys ) {
74 18 : AdjacencyMatrixBase::registerKeywords( keys );
75 36 : keys.add("atoms","BRIDGING_ATOMS","The list of atoms that can form the bridge between the two interesting parts "
76 : "of the structure.");
77 36 : keys.add("optional","SWITCH","The parameters of the two switchingfunction in the above formula");
78 36 : keys.add("optional","SWITCHA","The switchingfunction on the distance between bridging atoms and the atoms in "
79 : "group A");
80 36 : keys.add("optional","SWITCHB","The switchingfunction on the distance between the bridging atoms and the atoms in "
81 : "group B");
82 18 : }
83 :
84 8 : BridgeMatrix::BridgeMatrix(const ActionOptions&ao):
85 : Action(ao),
86 8 : AdjacencyMatrixBase(ao) {
87 : bool oneswitch;
88 : std::string sfinput,errors;
89 16 : parse("SWITCH",sfinput);
90 8 : if( sfinput.length()>0 ) {
91 8 : sf1.set(sfinput,errors);
92 8 : oneswitch=true;
93 8 : if( errors.length()!=0 ) {
94 0 : error("problem reading SWITCH keyword : " + errors );
95 : }
96 8 : sf2.set(sfinput,errors);
97 8 : if( errors.length()!=0 ) {
98 0 : error("problem reading SWITCH keyword : " + errors );
99 : }
100 : } else {
101 0 : parse("SWITCHA",sfinput);
102 0 : if(sfinput.length()>0) {
103 0 : sf1.set(sfinput,errors);
104 0 : oneswitch=false;
105 0 : if( errors.length()!=0 ) {
106 0 : error("problem reading SWITCHA keyword : " + errors );
107 : }
108 : sfinput.clear();
109 0 : parse("SWITCHB",sfinput);
110 0 : if(sfinput.length()==0) {
111 0 : error("found SWITCHA keyword without SWITCHB");
112 : }
113 0 : sf2.set(sfinput,errors);
114 0 : if( errors.length()!=0 ) {
115 0 : error("problem reading SWITCHB keyword : " + errors );
116 : }
117 : } else {
118 0 : error("missing definition of switching functions");
119 : }
120 : }
121 8 : log.printf(" distance between bridging atoms and atoms in GROUPA must be less than %s\n",sf1.description().c_str());
122 8 : log.printf(" distance between bridging atoms and atoms in GROUPB must be less than %s\n",sf2.description().c_str());
123 :
124 : // Setup link cells
125 8 : setLinkCellCutoff( oneswitch, sf1.get_dmax() + sf2.get_dmax() );
126 :
127 : // And check everything has been read in correctly
128 8 : checkRead();
129 8 : }
130 :
131 537 : double BridgeMatrix::calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const {
132 : double tot=0;
133 537 : if( pos2.modulo2()<epsilon ) {
134 : return 0.0;
135 : }
136 21359 : for(unsigned i=0; i<natoms; ++i) {
137 20822 : Vector dij= getPosition(i,myvals);
138 20822 : double dijm = dij.modulo2();
139 20822 : double dw1, w1=sf1.calculateSqr( dijm, dw1 );
140 20822 : if( dijm<epsilon ) {
141 : w1=0.0;
142 0 : dw1=0.0;
143 : }
144 20822 : Vector dik=pbcDistance( getPosition(i,myvals), pos2 );
145 20822 : double dikm=dik.modulo2();
146 20822 : double dw2, w2=sf2.calculateSqr( dikm, dw2 );
147 20822 : if( dikm<epsilon ) {
148 : w2=0.0;
149 0 : dw2=0.0;
150 : }
151 :
152 20822 : tot += w1*w2;
153 : // And finish the calculation
154 20822 : addAtomDerivatives( 0, -w2*dw1*dij, myvals );
155 20822 : addAtomDerivatives( 1, w1*dw2*dik, myvals );
156 20822 : addThirdAtomDerivatives( i, -w1*dw2*dik+w2*dw1*dij, myvals );
157 20822 : addBoxDerivatives( w1*(-dw2)*Tensor(dik,dik)+w2*(-dw1)*Tensor(dij,dij), myvals );
158 : }
159 537 : return tot;
160 : }
161 :
162 : }
163 : }
|
