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