Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 "DistanceFromContourBase.h"
23 :
24 : namespace PLMD {
25 : namespace contour {
26 :
27 12 : void DistanceFromContourBase::registerKeywords( Keywords& keys ) {
28 12 : Action::registerKeywords( keys );
29 12 : ActionWithValue::registerKeywords( keys );
30 12 : ActionAtomistic::registerKeywords( keys );
31 12 : ActionWithArguments::registerKeywords( keys );
32 12 : keys.remove("NUMERICAL_DERIVATIVES");
33 24 : keys.addInputKeyword("optional","ARG","vector","the label of the weights to use when constructing the density. If this keyword is not here the weights are assumed to be one.");
34 12 : keys.add("atoms","POSITIONS","the positions of the atoms that we are calculating the contour from");
35 12 : keys.add("atoms","ATOM","The atom whose perpendicular distance we are calculating from the contour");
36 12 : keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
37 12 : keys.add("compulsory","KERNEL","GAUSSIAN","the kernel function you are using. More details on the kernels available "
38 : "in plumed can be found in the documentation for KDE.");
39 12 : keys.add("compulsory","CUTOFF","6.25","the cutoff at which to stop evaluating the kernel functions is set equal to sqrt(2*x)*bandwidth in each direction where x is this number");
40 12 : keys.add("compulsory","CONTOUR","the value we would like for the contour");
41 12 : }
42 :
43 4 : DistanceFromContourBase::DistanceFromContourBase( const ActionOptions& ao ):
44 : Action(ao),
45 : ActionWithValue(ao),
46 : ActionAtomistic(ao),
47 : ActionWithArguments(ao),
48 : mymin(this),
49 4 : nactive(0) {
50 4 : if( getNumberOfArguments()>1 ) {
51 0 : error("should only use one argument for this action");
52 : }
53 4 : if( getNumberOfArguments()==1 && getPntrToArgument(0)->getRank()!=1 ) {
54 0 : error("ARG for distance from contour should be rank one");
55 : }
56 :
57 : // Read in the multicolvar/atoms
58 4 : std::vector<AtomNumber> atoms;
59 8 : parseAtomList("POSITIONS",atoms);
60 4 : if( atoms.empty() ) {
61 1 : error("no atoms were specified with POSITIONS");
62 : }
63 3 : std::vector<AtomNumber> origin;
64 6 : parseAtomList("ATOM",origin);
65 3 : if( origin.size()!=1 ) {
66 0 : error("should only specify one atom for origin keyword");
67 : }
68 3 : std::vector<AtomNumber> center;
69 3 : if( keywords.exists("ORIGIN") ) {
70 4 : parseAtomList("ORIGIN",center);
71 2 : if( center.size()!=1 ) {
72 0 : error("should only specify one atom for center keyword");
73 : }
74 : }
75 :
76 3 : if( center.size()==1 ) {
77 2 : log.printf(" calculating distance between atom %d and contour along vector connecting it to atom %d \n", origin[0].serial(),center[0].serial() );
78 : } else {
79 1 : log.printf(" calculating distance between atom %d and contour \n", origin[0].serial() );
80 : }
81 :
82 3 : log.printf(" contour is in field constructed from positions of atoms : ");
83 1028 : for(unsigned i=0; i<atoms.size(); ++i) {
84 1025 : log.printf("%d ",atoms[i].serial() );
85 : }
86 3 : if( getNumberOfArguments()==1 ) {
87 1 : if( getPntrToArgument(0)->getShape()[0]!=atoms.size() ) {
88 0 : error("mismatch between number of atoms and size of vector specified using ARG keyword");
89 : }
90 1 : log.printf("\n and weights from %s \n", getPntrToArgument(0)->getName().c_str() );
91 : } else {
92 2 : log.printf("\n all weights are set equal to one \n");
93 : }
94 : // Request everything we need
95 3 : active_list.resize( atoms.size(), 0 );
96 3 : std::vector<Value*> args( getArguments() );
97 3 : atoms.push_back( origin[0] );
98 3 : if( center.size()==1 ) {
99 2 : atoms.push_back( center[0] );
100 : }
101 3 : requestArguments( args );
102 3 : requestAtoms( atoms );
103 : // Fix to request arguments
104 3 : if( args.size()==1 ) {
105 1 : addDependency( args[0]->getPntrToAction() );
106 : }
107 :
108 : // Read in details of phase field construction
109 3 : parseVector("BANDWIDTH",bw);
110 3 : parse("KERNEL",kerneltype);
111 6 : parse("CONTOUR",contour);
112 : std::string errors;
113 6 : switchingFunction.set( kerneltype + " R_0=1.0 NOSTRETCH", errors );
114 3 : if( errors.length()!=0 ) {
115 0 : error("problem reading switching function description " + errors);
116 : }
117 : double det=1;
118 12 : for(unsigned i=0; i<bw.size(); ++i) {
119 9 : det*=bw[i]*bw[i];
120 : }
121 3 : gvol=1.0;
122 3 : if( kerneltype=="GAUSSIAN" ) {
123 3 : gvol=pow( 2*pi, 0.5*bw.size() ) * pow( det, 0.5 );
124 : }
125 3 : log.printf(" constructing phase field using %s kernels with bandwidth (%f, %f, %f) \n",kerneltype.c_str(), bw[0], bw[1], bw[2] );
126 :
127 : // And a cutoff
128 3 : std::vector<double> pp( bw.size(),0 );
129 : double dp2cutoff;
130 3 : parse("CUTOFF",dp2cutoff);
131 3 : double rcut = sqrt(2*dp2cutoff)*bw[0];
132 9 : for(unsigned j=1; j<bw.size(); ++j) {
133 6 : if( sqrt(2*dp2cutoff)*bw[j]>rcut ) {
134 0 : rcut=sqrt(2*dp2cutoff)*bw[j];
135 : }
136 : }
137 3 : rcut2=rcut*rcut;
138 :
139 : // Create the vector of values that holds the position
140 3 : forcesToApply.resize( 3*getNumberOfAtoms() + 9 );
141 15 : }
142 :
143 171 : void DistanceFromContourBase::lockRequests() {
144 : ActionWithArguments::lockRequests();
145 : ActionAtomistic::lockRequests();
146 171 : }
147 :
148 171 : void DistanceFromContourBase::unlockRequests() {
149 : ActionWithArguments::unlockRequests();
150 : ActionAtomistic::unlockRequests();
151 171 : }
152 :
153 267305 : double DistanceFromContourBase::evaluateKernel( const Vector& cpos, const Vector& apos, std::vector<double>& der ) const {
154 267305 : Vector distance = pbcDistance( getPosition(getNumberOfAtoms()-1), cpos );
155 : Vector dist2 = pbcDistance( distance, apos );
156 : double dval=0;
157 1069220 : for(unsigned j=0; j<3; ++j) {
158 801915 : der[j] = dist2[j]/bw[j];
159 801915 : dval += der[j]*der[j];
160 : }
161 267305 : double dfunc, newval = switchingFunction.calculateSqr( dval, dfunc ) / gvol;
162 267305 : double tmp = dfunc / gvol;
163 1069220 : for(unsigned j=0; j<3; ++j) {
164 801915 : der[j] *= -tmp;
165 : }
166 267305 : return newval;
167 : }
168 :
169 3458 : double DistanceFromContourBase::getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der ) {
170 : // Transer the position to the local Vector
171 13832 : for(unsigned j=0; j<3; ++j) {
172 10374 : pval[j] = x[j];
173 : }
174 : // Now find the contour
175 : double sumk = 0, sumd = 0;
176 3458 : std::vector<double> pp(3), ddd(3,0);
177 270489 : for(unsigned i=0; i<nactive; ++i) {
178 267031 : double newval = evaluateKernel( getPosition(active_list[i]), pval, ddd );
179 :
180 267031 : if( getNumberOfArguments()==1 ) {
181 174590 : sumk += getPntrToArgument(0)->get(active_list[i])*newval;
182 174590 : sumd += newval;
183 : } else {
184 92441 : sumk += newval;
185 : }
186 : }
187 3458 : if( getNumberOfArguments()==0 ) {
188 3016 : return sumk - contour;
189 : }
190 442 : if( fabs(sumk)<epsilon ) {
191 207 : return -contour;
192 : }
193 235 : return (sumk/sumd) - contour;
194 3458 : }
195 :
196 171 : void DistanceFromContourBase::apply() {
197 171 : if( !checkForForces() ) {
198 34 : return ;
199 : }
200 137 : unsigned ind=0;
201 137 : setForcesOnAtoms( getForcesToApply(), ind );
202 : }
203 :
204 : }
205 : }
|
