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 "Gradient.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/HistogramBead.h"
26 :
27 : namespace PLMD {
28 : namespace crystallization {
29 :
30 : //+PLUMEDOC MCOLVARF GRADIENT
31 : /*
32 : Calculate the gradient of the average value of a multicolvar value
33 :
34 : This command allows you to calculate the collective variable discussed in \cite fede-grad.
35 :
36 : \par Examples
37 :
38 : The input below calculates the gradient of the density of atoms in the manner
39 : described in \cite fede-grad in order to detect whether or not atoms are distributed
40 : uniformly along the x-axis of the simulation cell.
41 :
42 : \plumedfile
43 : d1: DENSITY SPECIES=1-50
44 : s1: GRADIENT ORIGIN=1 DATA=d1 DIR=x NBINS=4 SIGMA=1.0
45 : PRINT ARG=s1 FILE=colvar
46 : \endplumedfile
47 :
48 : The input below calculates the coordination numbers of the 50 atoms in the simulation cell.
49 : The gradient of this quantity is then evaluated in the manner described using the equation above
50 : to detect whether the average values of the coordination number are uniformly distributed along the
51 : x-axis of the simulation cell.
52 :
53 : \plumedfile
54 : d2: COORDINATIONNUMBER SPECIES=1-50 SWITCH={RATIONAL R_0=2.0} MORE_THAN={EXP R_0=4.0}
55 : s2: GRADIENT ORIGIN=1 DATA=d2 DIR=x NBINS=4 SIGMA=1.0
56 : PRINT ARG=s2 FILE=colvar
57 : \endplumedfile
58 :
59 : */
60 : //+ENDPLUMEDOC
61 :
62 13789 : PLUMED_REGISTER_ACTION(Gradient,"GRADIENT")
63 :
64 6 : void Gradient::registerKeywords( Keywords& keys ) {
65 6 : VolumeGradientBase::registerKeywords( keys );
66 12 : keys.add("atoms","ORIGIN","we will use the position of this atom as the origin in our calculation");
67 12 : keys.add("compulsory","DIR","xyz","the directions in which we are calculating the gradient. Should be x, y, z, xy, xz, yz or xyz");
68 12 : keys.add("compulsory","NBINS","number of bins to use in each direction for the calculation of the gradient");
69 12 : keys.add("compulsory","SIGMA","1.0","the width of the function to be used for kernel density estimation");
70 12 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
71 6 : }
72 :
73 2 : Gradient::Gradient(const ActionOptions&ao):
74 : Action(ao),
75 : VolumeGradientBase(ao),
76 2 : nbins(3) {
77 : std::vector<AtomNumber> atom;
78 4 : parseAtomList("ORIGIN",atom);
79 2 : if( atom.size()!=1 ) {
80 0 : error("should only be one atom specified");
81 : }
82 2 : log.printf(" origin is at position of atom : %d\n",atom[0].serial() );
83 :
84 : std::string direction;
85 4 : parse("DIR",direction);
86 : std::vector<unsigned> tbins;
87 2 : parseVector("NBINS",tbins);
88 4 : for(unsigned i=0; i<tbins.size(); ++i) {
89 2 : if( tbins[i]<2 ) {
90 0 : error("Number of grid points should be greater than 1");
91 : }
92 : }
93 :
94 2 : if( direction=="x" ) {
95 2 : if( tbins.size()!=1 ) {
96 0 : error("mismatch between number of bins and direction");
97 : }
98 2 : nbins[0]=tbins[0];
99 2 : nbins[1]=0;
100 2 : nbins[2]=0;
101 0 : } else if( direction=="y" ) {
102 0 : if( tbins.size()!=1 ) {
103 0 : error("mismatch between number of bins and direction");
104 : }
105 0 : nbins[0]=0;
106 0 : nbins[1]=tbins[0];
107 0 : nbins[2]=0;
108 0 : } else if( direction=="z" ) {
109 0 : if( tbins.size()!=1 ) {
110 0 : error("mismatch between number of bins and direction");
111 : }
112 0 : nbins[0]=0;
113 0 : nbins[1]=0;
114 0 : nbins[2]=tbins[0];
115 0 : } else if( direction=="xy" ) {
116 0 : if( tbins.size()!=2 ) {
117 0 : error("mismatch between number of bins and direction");
118 : }
119 0 : nbins[0]=tbins[0];
120 0 : nbins[1]=tbins[1];
121 0 : nbins[2]=0;
122 0 : } else if( direction=="xz" ) {
123 0 : if( tbins.size()!=2 ) {
124 0 : error("mismatch between number of bins and direction");
125 : }
126 0 : nbins[0]=tbins[0];
127 0 : nbins[1]=0;
128 0 : nbins[2]=tbins[1];
129 0 : } else if( direction=="yz" ) {
130 0 : if( tbins.size()!=2 ) {
131 0 : error("mismatch between number of bins and direction");
132 : }
133 0 : nbins[0]=0;
134 0 : nbins[1]=tbins[0];
135 0 : nbins[2]=tbins[1];
136 0 : } else if( direction=="xyz" ) {
137 0 : if( tbins.size()!=3 ) {
138 0 : error("mismatch between number of bins and direction");
139 : }
140 0 : nbins[0]=tbins[0];
141 0 : nbins[1]=tbins[1];
142 0 : nbins[2]=tbins[2];
143 : } else {
144 0 : error( direction + " is not valid gradient direction");
145 : }
146 :
147 : // Find number of quantities
148 2 : if( getPntrToMultiColvar()->isDensity() ) {
149 1 : vend=2;
150 1 : } else if( getPntrToMultiColvar()->getNumberOfQuantities()==2 ) {
151 1 : vend=2;
152 : } else {
153 0 : vend = getPntrToMultiColvar()->getNumberOfQuantities();
154 : }
155 2 : nquantities = vend + nbins[0] + nbins[1] + nbins[2];
156 :
157 : // Output some nice information
158 2 : std::string functype=getPntrToMultiColvar()->getName();
159 2 : std::transform( functype.begin(), functype.end(), functype.begin(), [](unsigned char c) {
160 25 : return std::tolower(c);
161 : } );
162 2 : log.printf(" calculating gradient of %s in %s direction \n",functype.c_str(), direction.c_str() );
163 4 : log<<" Bibliography:"<<plumed.cite("Giberti, Tribello and Parrinello, J. Chem. Theory Comput., 9, 2526 (2013)")<<"\n";
164 :
165 2 : parse("SIGMA",sigma);
166 2 : parse("KERNEL",kerneltype);
167 2 : checkRead();
168 2 : requestAtoms(atom);
169 :
170 : // And setup the vessel
171 : std::string input;
172 2 : addVessel( "GRADIENT", input, -1 );
173 : // And resize everything
174 2 : readVesselKeywords();
175 2 : }
176 :
177 232 : void Gradient::setupRegions() {
178 : // if( !getPbc().isOrthorombic() ) error("cell must be orthorhombic when using gradient");
179 232 : }
180 :
181 11600 : void Gradient::calculateAllVolumes( const unsigned& curr, MultiValue& outvals ) const {
182 : // Setup the bead
183 11600 : HistogramBead bead;
184 : bead.isNotPeriodic();
185 11600 : bead.setKernelType( kerneltype );
186 :
187 11600 : Vector cpos = pbcDistance( getPosition(0), getPntrToMultiColvar()->getCentralAtomPos( curr ) );
188 : // Note we use the pbc from base multicolvar so that we get numerical derivatives correct
189 11600 : Vector oderiv, fpos = (getPntrToMultiColvar()->getPbc()).realToScaled( cpos );
190 :
191 11600 : Vector deriv;
192 11600 : unsigned nbase=vend;
193 11600 : std::vector<Vector> refder(1);
194 11600 : Tensor vir;
195 11600 : vir.zero();
196 46400 : for(unsigned idir=0; idir<3; ++idir) {
197 34800 : deriv[0]=deriv[1]=deriv[2]=0.0;
198 34800 : double delx = 1.0 / static_cast<double>( nbins[idir] );
199 81200 : for(unsigned jbead=0; jbead<nbins[idir]; ++jbead) {
200 : // Calculate what box we are in
201 46400 : bead.set( -0.5+jbead*delx, -0.5+(jbead+1)*delx, sigma );
202 46400 : double weight=bead.calculate( fpos[0], deriv[idir] );
203 46400 : oderiv = (getPntrToMultiColvar()->getPbc()).realToScaled( deriv );
204 : // Set and derivatives
205 46400 : refder[0]=-oderiv; // vir = -Tensor(cpos,oderiv);
206 46400 : setNumberInVolume( nbase+jbead, curr, weight, oderiv, vir, refder, outvals );
207 : // addReferenceAtomDerivatives( nbase+jbead, 0, -oderiv );
208 : // addBoxDerivatives( nbase+jbead, -Tensor(cpos,oderiv) );
209 : }
210 34800 : nbase += nbins[idir];
211 : }
212 11600 : }
213 :
214 : }
215 : }
|
