Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2018 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 : \verbatim
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 : \endverbatim
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 : \verbatim
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 : \endverbatim
58 :
59 : */
60 : //+ENDPLUMEDOC
61 :
62 2525 : PLUMED_REGISTER_ACTION(Gradient,"GRADIENT")
63 :
64 3 : void Gradient::registerKeywords( Keywords& keys ) {
65 3 : VolumeGradientBase::registerKeywords( keys );
66 3 : keys.add("atoms","ORIGIN","we will use the position of this atom as the origin in our calculation");
67 3 : keys.add("compulsory","DIR","xyz","the directions in which we are calculating the graident. Should be x, y, z, xy, xz, yz or xyz");
68 3 : keys.add("compulsory","NBINS","number of bins to use in each direction for the calculation of the gradient");
69 3 : keys.add("compulsory","SIGMA","1.0","the width of the function to be used for kernel density estimation");
70 3 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
71 3 : }
72 :
73 2 : Gradient::Gradient(const ActionOptions&ao):
74 : Action(ao),
75 : VolumeGradientBase(ao),
76 2 : nbins(3)
77 : {
78 2 : std::vector<AtomNumber> atom;
79 2 : parseAtomList("ORIGIN",atom);
80 2 : if( atom.size()!=1 ) error("should only be one atom specified");
81 2 : log.printf(" origin is at position of atom : %d\n",atom[0].serial() );
82 :
83 4 : std::string direction; parse("DIR",direction);
84 4 : std::vector<unsigned> tbins; parseVector("NBINS",tbins);
85 4 : for(unsigned i=0; i<tbins.size(); ++i) {
86 2 : if( tbins[i]<2 ) error("Number of grid points should be greater than 1");
87 : }
88 :
89 2 : if( direction=="x" ) {
90 2 : if( tbins.size()!=1 ) error("mismatch between number of bins and direction");
91 2 : nbins[0]=tbins[0]; nbins[1]=0; nbins[2]=0;
92 0 : } else if( direction=="y" ) {
93 0 : if( tbins.size()!=1 ) error("mismatch between number of bins and direction");
94 0 : nbins[0]=0; nbins[1]=tbins[0]; nbins[2]=0;
95 0 : } else if( direction=="z" ) {
96 0 : if( tbins.size()!=1 ) error("mismatch between number of bins and direction");
97 0 : nbins[0]=0; nbins[1]=0; nbins[2]=tbins[0];
98 0 : } else if( direction=="xy" ) {
99 0 : if( tbins.size()!=2 ) error("mismatch between number of bins and direction");
100 0 : nbins[0]=tbins[0]; nbins[1]=tbins[1]; nbins[2]=0;
101 0 : } else if( direction=="xz" ) {
102 0 : if( tbins.size()!=2 ) error("mismatch between number of bins and direction");
103 0 : nbins[0]=tbins[0]; nbins[1]=0; nbins[2]=tbins[1];
104 0 : } else if( direction=="yz" ) {
105 0 : if( tbins.size()!=2 ) error("mismatch between number of bins and direction");
106 0 : nbins[0]=0; nbins[1]=tbins[0]; nbins[2]=tbins[1];
107 0 : } else if( direction=="xyz" ) {
108 0 : if( tbins.size()!=3 ) error("mismatch between number of bins and direction");
109 0 : nbins[0]=tbins[0]; nbins[1]=tbins[1]; nbins[2]=tbins[2];
110 : } else {
111 0 : error( direction + " is not valid gradient direction");
112 : }
113 :
114 : // Find number of quantities
115 2 : if( getPntrToMultiColvar()->isDensity() ) vend=2;
116 1 : else if( getPntrToMultiColvar()->getNumberOfQuantities()==2 ) vend=2;
117 0 : else vend = getPntrToMultiColvar()->getNumberOfQuantities();
118 2 : nquantities = vend + nbins[0] + nbins[1] + nbins[2];
119 :
120 : // Output some nice information
121 4 : std::string functype=getPntrToMultiColvar()->getName();
122 2 : std::transform( functype.begin(), functype.end(), functype.begin(), tolower );
123 2 : log.printf(" calculating gradient of %s in %s direction \n",functype.c_str(), direction.c_str() );
124 2 : log<<" Bibliography:"<<plumed.cite("Giberti, Tribello and Parrinello, J. Chem. Theory Comput., 9, 2526 (2013)")<<"\n";
125 :
126 2 : parse("SIGMA",sigma); parse("KERNEL",kerneltype);
127 2 : checkRead(); requestAtoms(atom);
128 :
129 : // And setup the vessel
130 4 : std::string input; addVessel( "GRADIENT", input, -1 );
131 : // And resize everything
132 4 : readVesselKeywords();
133 2 : }
134 :
135 232 : void Gradient::setupRegions() {
136 : // if( !getPbc().isOrthorombic() ) error("cell must be orthorhombic when using gradient");
137 232 : }
138 :
139 11530 : void Gradient::calculateAllVolumes( const unsigned& curr, MultiValue& outvals ) const {
140 : // Setup the bead
141 11530 : HistogramBead bead; bead.isNotPeriodic(); bead.setKernelType( kerneltype );
142 :
143 11574 : Vector cpos = pbcDistance( getPosition(0), getPntrToMultiColvar()->getCentralAtomPos( curr ) );
144 : // Note we use the pbc from base multicolvar so that we get numerical derivatives correct
145 11516 : Vector oderiv, fpos = (getPntrToMultiColvar()->getPbc()).realToScaled( cpos );
146 :
147 11587 : Vector deriv; unsigned nbase=vend; std::vector<Vector> refder(1); Tensor vir; vir.zero();
148 46348 : for(unsigned idir=0; idir<3; ++idir) {
149 34750 : deriv[0]=deriv[1]=deriv[2]=0.0;
150 34764 : double delx = 1.0 / static_cast<double>( nbins[idir] );
151 81044 : for(unsigned jbead=0; jbead<nbins[idir]; ++jbead) {
152 : // Calculate what box we are in
153 46346 : bead.set( -0.5+jbead*delx, -0.5+(jbead+1)*delx, sigma );
154 46279 : double weight=bead.calculate( fpos[0], deriv[idir] );
155 46362 : oderiv = (getPntrToMultiColvar()->getPbc()).realToScaled( deriv );
156 : // Set and derivatives
157 46391 : refder[0]=-oderiv; // vir = -Tensor(cpos,oderiv);
158 46387 : setNumberInVolume( nbase+jbead, curr, weight, oderiv, vir, refder, outvals );
159 : // addReferenceAtomDerivatives( nbase+jbead, 0, -oderiv );
160 : // addBoxDerivatives( nbase+jbead, -Tensor(cpos,oderiv) );
161 : }
162 34765 : nbase += nbins[idir];
163 11598 : }
164 11600 : }
165 :
166 : }
167 2523 : }
|
