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