Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2020 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 : #ifndef __PLUMED_volumes_VolumeInCylinder_h
23 : #define __PLUMED_volumes_VolumeInCylinder_h
24 : #include "tools/Pbc.h"
25 : #include "tools/SwitchingFunction.h"
26 : #include "ActionVolume.h"
27 : #include "tools/HistogramBead.h"
28 :
29 :
30 : namespace PLMD {
31 : namespace volumes {
32 :
33 3 : class VolumeInCylinder {
34 : public:
35 : using precision=double;
36 : typedef VolumeIn<precision> VolumeInput;
37 : typedef VolumeOut<precision> VolumeOutput;
38 : bool docylinder;
39 : double min, max, sigma;
40 : HistogramBead::KernelType kerneltype;
41 : std::array<unsigned,3> dir;
42 : #ifdef __PLUMED_HAS_OPENACC
43 : SwitchingFunctionAccelerable switchingFunction;
44 : #else
45 : SwitchingFunction switchingFunction;
46 : #endif //__PLUMED_HAS_OPENACC
47 : static void registerKeywords( Keywords& keys );
48 : template<typename TPM>
49 : void parseInput( ActionVolume<VolumeInCylinder, TPM>* action );
50 : template<typename TPM>
51 : void setupRegions( ActionVolume<VolumeInCylinder, TPM>* action,
52 : const Pbc& pbc,
53 : const std::vector<Vector>& positions ) {}
54 : template<typename TPM>
55 : static void parseAtoms( ActionVolume<VolumeInCylinder, TPM>* action,
56 : std::vector<AtomNumber>& atom );
57 : static void calculateNumberInside( const VolumeInput& input,
58 : const VolumeInCylinder& actioninput,
59 : VolumeOutput& output );
60 : #ifdef __PLUMED_HAS_OPENACC
61 : void toACCDevice() const {
62 : #pragma acc enter data copyin(this[0:1], \
63 : docylinder, min, max, sigma, kerneltype, dir[0:3])
64 : switchingFunction.toACCDevice();
65 : }
66 : void removeFromACCDevice() const {
67 : switchingFunction.removeFromACCDevice();
68 : #pragma acc exit data delete(dir[0:3], kerneltype, sigma, max, min, \
69 : docylinder, this[0:1])
70 : }
71 : #endif //__PLUMED_HAS_OPENACC
72 : };
73 :
74 11 : void VolumeInCylinder::registerKeywords( Keywords& keys ) {
75 22 : keys.setDisplayName("INCYLINDER");
76 11 : keys.add("atoms","CENTER","the atom whose vicinity we are interested in examining");
77 11 : keys.add("optional","SIGMA","the width of the function to be used for kernel density estimation");
78 11 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
79 11 : keys.add("compulsory","DIRECTION","the direction of the long axis of the cylinder. Must be x, y or z");
80 11 : keys.add("compulsory","RADIUS","a switching function that gives the extent of the cylinder in the plane perpendicular to the direction");
81 11 : keys.add("compulsory","LOWER","0.0","the lower boundary on the direction parallel to the long axis of the cylinder");
82 11 : keys.add("compulsory","UPPER","0.0","the upper boundary on the direction parallel to the long axis of the cylinder");
83 22 : keys.linkActionInDocs("RADIUS","LESS_THAN");
84 11 : }
85 :
86 : template<typename TPM>
87 1 : void VolumeInCylinder::parseInput( ActionVolume<VolumeInCylinder, TPM>* action ) {
88 2 : action->parse("SIGMA",sigma);
89 : std::string mykerneltype;
90 1 : action->parse("KERNEL",mykerneltype);
91 1 : kerneltype=HistogramBead::getKernelType(mykerneltype);
92 : std::string sdir;
93 1 : action->parse("DIRECTION",sdir);
94 1 : if( sdir=="X") {
95 0 : dir[0]=1;
96 0 : dir[1]=2;
97 0 : dir[2]=0;
98 1 : } else if( sdir=="Y") {
99 0 : dir[0]=0;
100 0 : dir[1]=2;
101 0 : dir[2]=1;
102 1 : } else if( sdir=="Z") {
103 1 : dir[0]=0;
104 1 : dir[1]=1;
105 1 : dir[2]=2;
106 : } else {
107 0 : action->error(sdir + "is not a valid direction. Should be X, Y or Z");
108 : }
109 1 : action->log.printf(" cylinder's long axis is along %s axis\n",sdir.c_str() );
110 :
111 : std::string errors;
112 : std::string swinput;
113 2 : action->parse("RADIUS",swinput);
114 1 : if(swinput.length()==0) {
115 0 : action->error("missing RADIUS keyword");
116 : }
117 1 : switchingFunction.set(swinput,errors);
118 1 : if( errors.length()!=0 ) {
119 0 : action->error("problem reading RADIUS keyword : " + errors );
120 : }
121 1 : action->log.printf(" radius of cylinder is given by %s \n", ( switchingFunction.description() ).c_str() );
122 :
123 1 : docylinder=false;
124 1 : action->parse("LOWER",min);
125 1 : action->parse("UPPER",max);
126 1 : if( min!=0.0 || max!=0.0 ) {
127 1 : if( min>max ) {
128 0 : action->error("minimum of cylinder should be less than maximum");
129 : }
130 1 : docylinder=true;
131 1 : action->log.printf(" cylinder extends from %f to %f along the %s axis\n",min,max,sdir.c_str() );
132 : }
133 1 : }
134 :
135 : template<typename TPM>
136 1 : void VolumeInCylinder::parseAtoms( ActionVolume<VolumeInCylinder, TPM>* action,
137 : std::vector<AtomNumber>& atom ) {
138 2 : action->parseAtomList("CENTER",atom);
139 1 : if( atom.size()!=1 ) {
140 0 : action->error("should only be one atom specified");
141 : }
142 1 : action->log.printf(" center of cylinder is at position of atom : %d\n",atom[0].serial() );
143 1 : }
144 :
145 8000 : void VolumeInCylinder::calculateNumberInside( const VolumeInput& input,
146 : const VolumeInCylinder& actioninput,
147 : VolumeOutput& output ) {
148 : // Calculate position of atom wrt to origin
149 8000 : Vector fpos=input.pbc.distance( Vector(input.refpos[0][0],input.refpos[0][1],input.refpos[0][2]), Vector(input.cpos[0],input.cpos[1],input.cpos[2]) );
150 :
151 : double vcylinder, dcylinder;
152 8000 : if( actioninput.docylinder ) {
153 8000 : HistogramBead bead( actioninput.kerneltype,
154 8000 : actioninput.min, actioninput.max, actioninput.sigma );
155 8000 : vcylinder=bead.calculate( fpos[actioninput.dir[2]], dcylinder );
156 : } else {
157 : vcylinder=1.0;
158 0 : dcylinder=0.0;
159 : }
160 :
161 8000 : const double dd = fpos[actioninput.dir[0]]*fpos[actioninput.dir[0]] + fpos[actioninput.dir[1]]*fpos[actioninput.dir[1]];
162 : double dfunc;
163 8000 : double vswitch = actioninput.switchingFunction.calculateSqr( dd, dfunc );
164 8000 : output.values[0]=vswitch*vcylinder;
165 8000 : output.derivatives[actioninput.dir[0]]=vcylinder*dfunc*fpos[actioninput.dir[0]];
166 8000 : output.derivatives[actioninput.dir[1]]=vcylinder*dfunc*fpos[actioninput.dir[1]];
167 8000 : output.derivatives[actioninput.dir[2]]=vswitch*dcylinder;
168 : // Add derivatives wrt to position of origin atom
169 8000 : output.refders[0][0] = -output.derivatives[0];
170 8000 : output.refders[0][1] = -output.derivatives[1];
171 8000 : output.refders[0][2] = -output.derivatives[2];
172 : // Add virial contribution
173 16000 : output.virial.set( 0, -Tensor(fpos,Vector(output.derivatives[0], output.derivatives[1], output.derivatives[2])) );
174 8000 : }
175 :
176 : }
177 : }
178 : #endif //__PLUMED_volumes_VolumeInCylinder_h
|
