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