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