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