Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 "OrientationSphere.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/Torsion.h"
25 : #include "tools/KernelFunctions.h"
26 :
27 : //+PLUMEDOC MCOLVARF SMAC
28 : /*
29 : Calculate a variant on the SMAC collective variable
30 :
31 : This variable is discussed in \cite smac-paper
32 :
33 : The SMAC collective variable can be used to study the formation of molecular solids
34 : from either the melt or from solution. The idea behind this variable is that what
35 : differentiates a molecular solid from a molecular liquid is an alignment of
36 : internal vectors in neighboring molecules. In other words, the relative orientation
37 : of neighboring molecules is no longer random as it is in a liquid. In a solid particular
38 : torsional angles between molecules are preferred. As such this CV calculates the following
39 : average:
40 :
41 : \f[
42 : s_i = \frac{ \left\{ 1 - \psi\left[ \sum_{j \ne i} \sigma(r_{ij}) \right] \right\} \sum_{j \ne i} \sigma(r_{ij}) \sum_n K_n(\theta_{ij}) }{ \sum_{j \ne i} \sigma(r_{ij}) }
43 : \f]
44 :
45 : In this expression \f$r_{ij}\f$ is the distance between molecule \f$i\f$ and molecule \f$j\f$ and \f$\sigma(r_{ij})\f$ is a
46 : \ref switchingfunction that acts on this distance. By including this switching function in the second summation in the
47 : numerator and in the denominator we are thus ensuring that we calculate an average over the molecules in the first coordination
48 : sphere of molecule \f$i\f$. All molecules in higher coordination sphere will essentially contribute zero to the sums in the
49 : above expression because their \f$\sigma(r_{ij})\f$ will be very small. \f$\psi\f$ is also a switching function. The term
50 : including \f$\psi\f$ in the numerator is there to ensure that only those molecules that are attached to a reasonably large
51 : number of molecules. It is important to include this "more than" switching function when you are simulating nucleation
52 : from solution with this CV. Lastly, the $K_n functions are \ref kernelfunctions that take the torsion angle, \f$\theta_{ij}\f$, between the
53 : internal orientation vectors for molecules \f$i\f$ and \f$j\f$ as input. These kernel functions should be set so that they are
54 : equal to one when the relative orientation of the molecules are as they are in the solid and equal to zero otherwise.
55 : The final \f$s_i\f$ quantity thus measures whether (on average) the molecules in the first coordination sphere around molecule \f$i\f$
56 : are oriented as they would be in the solid. Furthermore, this Action is a multicolvar so you can calculate the \f$s_i\f$ values
57 : for all the molecules in your system simultaneously and then determine the average, the number less than and so on.
58 :
59 : \par Examples
60 :
61 : In the example below the orientation of the molecules in the system is determined by calculating the
62 : vector that connects a pair of atoms. SMAC is then used to determine whether the molecules are sitting
63 : in a solid or liquid like environment. We can determine whether the environment is solid or liquid like because in the solid the torsional angle between
64 : the bond vectors on adjacent molecules is close to 0 or \f$\pi\f$. The final quantity that is output to the colvar
65 : file measures the number of molecules that have a SMAC parameter that is greater than 0.7. N.B. By using
66 : the indices of three atoms for each of the MOL keywords below we are telling PLUMED to use the first two
67 : numbers to determine the orientation of the molecule that will ultimately be used when calculating the \f$\theta_{ij}\f$
68 : terms in the formula above. The atom with the third index meanwhile is used when we calculate \f$r_{ij}\f$.
69 :
70 : \plumedfile
71 : MOLECULES ...
72 : MOL1=9,10,9
73 : MOL2=89,90,89
74 : MOL3=473,474,473
75 : MOL4=1161,1162,1161
76 : MOL5=1521,1522,1521
77 : MOL6=1593,1594,1593
78 : MOL7=1601,1602,1601
79 : MOL8=2201,2202,2201
80 : LABEL=m3
81 : ... MOLECULES
82 :
83 : SMAC ...
84 : SPECIES=m3 LOWMEM
85 : KERNEL1={GAUSSIAN CENTER=0 SIGMA=0.480} KERNEL2={GAUSSIAN CENTER=pi SIGMA=0.480}
86 : SWITCH={RATIONAL R_0=0.6} MORE_THAN={RATIONAL R_0=0.7} SWITCH_COORD={EXP R_0=4}
87 : LABEL=s2
88 : ... SMAC
89 :
90 : PRINT ARG=s2.* FILE=colvar
91 : \endplumedfile
92 :
93 : This second example works in a way that is very similar to the previous command. Now, however,
94 : the orientation of the molecules is determined by finding the plane that contains the positions
95 : of three atoms.
96 :
97 : \plumedfile
98 : PLANES ...
99 : MOL1=9,10,11
100 : MOL2=89,90,91
101 : MOL3=473,474,475
102 : MOL4=1161,1162,1163
103 : MOL5=1521,1522,1523
104 : MOL6=1593,1594,1595
105 : MOL7=1601,1602,1603
106 : MOL8=2201,2202,2203
107 : VMEAN
108 : LABEL=m3
109 : ... PLANES
110 :
111 : SMAC ...
112 : SPECIES=m3 LOWMEM
113 : KERNEL1={GAUSSIAN CENTER=0 SIGMA=0.480} KERNEL2={GAUSSIAN CENTER=pi SIGMA=0.480}
114 : SWITCH={RATIONAL R_0=0.6} MORE_THAN={RATIONAL R_0=0.7} SWITCH_COORD={EXP R_0=3.0}
115 : LABEL=s2
116 : ... SMAC
117 :
118 : PRINT ARG=s2.* FILE=colvar
119 : \endplumedfile
120 :
121 : */
122 : //+ENDPLUMEDOC
123 :
124 : namespace PLMD {
125 : namespace crystallization {
126 :
127 : class SMAC : public OrientationSphere {
128 : private:
129 : std::vector<KernelFunctions> kernels;
130 : SwitchingFunction coord_switch;
131 : public:
132 : static void registerKeywords( Keywords& keys );
133 : explicit SMAC(const ActionOptions& ao);
134 : double computeVectorFunction( const Vector& conn, const std::vector<double>& vec1, const std::vector<double>& vec2,
135 : Vector& dconn, std::vector<double>& dvec1, std::vector<double>& dvec2 ) const override;
136 : double calculateCoordinationPrefactor( const double& coord, double& df ) const override;
137 : };
138 :
139 13795 : PLUMED_REGISTER_ACTION(SMAC,"SMAC")
140 :
141 9 : void SMAC::registerKeywords( Keywords& keys ) {
142 9 : OrientationSphere::registerKeywords(keys);
143 18 : keys.add("numbered","KERNEL","The kernels used in the function of the angle");
144 18 : keys.add("compulsory","SWITCH_COORD","This keyword is used to define the coordination switching function.");
145 18 : keys.reset_style("KERNEL","compulsory");
146 9 : }
147 :
148 5 : SMAC::SMAC(const ActionOptions& ao):
149 : Action(ao),
150 5 : OrientationSphere(ao) {
151 5 : if( mybasemulticolvars.size()==0 ) {
152 0 : error("SMAC must take multicolvar as input");
153 : }
154 11 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
155 6 : if( (mybasemulticolvars[i]->getNumberOfQuantities()-2)%3!=0 ) {
156 0 : error("SMAC is only possible with three dimensional vectors");
157 : }
158 : }
159 :
160 : std::string kernelinpt;
161 10 : for(int i=1;; i++) {
162 30 : if( !parseNumbered("KERNEL",i,kernelinpt) ) {
163 : break;
164 : }
165 10 : KernelFunctions mykernel( kernelinpt );
166 10 : kernels.push_back( mykernel );
167 10 : }
168 5 : if( kernels.size()==0 ) {
169 0 : error("no kernels defined");
170 : }
171 :
172 : std::string sw, errors;
173 10 : parse("SWITCH_COORD",sw);
174 5 : if(sw.length()==0) {
175 0 : error("SWITCH_COORD keyword is missing");
176 : }
177 5 : coord_switch.set(sw,errors);
178 5 : if(errors.length()>0) {
179 0 : error("the following errors were found in input to SWITCH_COORD : " + errors );
180 : }
181 :
182 5 : }
183 :
184 1025856 : double SMAC::computeVectorFunction( const Vector& conn, const std::vector<double>& vec1, const std::vector<double>& vec2,
185 : Vector& dconn, std::vector<double>& dvec1, std::vector<double>& dvec2 ) const {
186 :
187 1025856 : unsigned nvectors = ( vec1.size() - 2 ) / 3;
188 1025856 : plumed_assert( (vec1.size()-2)%3==0 );
189 1025856 : std::vector<Vector> dv1(nvectors), dv2(nvectors), tdconn(nvectors);
190 : Torsion t;
191 1025856 : std::vector<Vector> v1(nvectors), v2(nvectors);
192 : std::vector<std::unique_ptr<Value>> pos;
193 2051712 : for(unsigned i=0; i<nvectors; ++i) {
194 1025856 : pos.emplace_back( Tools::make_unique<Value>() );
195 2051712 : pos[i]->setDomain( "-pi", "pi" );
196 : }
197 :
198 2051712 : for(unsigned j=0; j<nvectors; ++j) {
199 4103424 : for(unsigned k=0; k<3; ++k) {
200 3077568 : v1[j][k]=vec1[2+3*j+k];
201 3077568 : v2[j][k]=vec2[2+3*j+k];
202 : }
203 1025856 : double angle = t.compute( v1[j], conn, v2[j], dv1[j], tdconn[j], dv2[j] );
204 : pos[j]->set( angle );
205 : }
206 :
207 1025856 : auto pos_ptr=Tools::unique2raw(pos);
208 :
209 : double ans=0;
210 1025856 : std::vector<double> deriv( nvectors ), df( nvectors, 0 );
211 3077568 : for(unsigned i=0; i<kernels.size(); ++i) {
212 2051712 : ans += kernels[i].evaluate( pos_ptr, deriv );
213 4103424 : for(unsigned j=0; j<nvectors; ++j) {
214 2051712 : df[j] += deriv[j];
215 : }
216 : }
217 1025856 : dconn.zero();
218 2051712 : for(unsigned j=0; j<nvectors; ++j) {
219 1025856 : dconn += df[j]*tdconn[j];
220 : }
221 2051712 : for(unsigned j=0; j<nvectors; ++j) {
222 4103424 : for(unsigned k=0; k<3; ++k) {
223 3077568 : dvec1[2+3*j+k]=df[j]*dv1[j][k];
224 3077568 : dvec2[2+3*j+k]=df[j]*dv2[j][k];
225 : }
226 : }
227 1025856 : return ans;
228 1025856 : }
229 :
230 4464 : double SMAC::calculateCoordinationPrefactor( const double& coord, double& df ) const {
231 4464 : double f=1-coord_switch.calculate( coord, df );
232 4464 : df*=-coord;
233 4464 : return f;
234 : }
235 :
236 : }
237 : }
|
