Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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_VolumeAround_h
23 : #define __PLUMED_volumes_VolumeAround_h
24 : #include "tools/Pbc.h"
25 : #include "tools/HistogramBead.h"
26 : #include "ActionVolume.h"
27 :
28 :
29 : namespace PLMD {
30 : namespace volumes {
31 :
32 98 : class VolumeAround {
33 : public:
34 : using precision=double;
35 : typedef VolumeIn<precision> VolumeInput;
36 : typedef VolumeOut<precision> VolumeOutput;
37 : bool dox{true}, doy{true}, doz{true};
38 : double sigma;
39 : double xlow{0.0}, xhigh{0.0};
40 : double ylow{0.0}, yhigh{0.0};
41 : double zlow{0.0}, zhigh{0.0};
42 : HistogramBead::KernelType kerneltype;
43 : static void registerKeywords( Keywords& keys );
44 : template<typename TPM>
45 : void parseInput( ActionVolume<VolumeAround, TPM>* action );
46 : template<typename TPM>
47 : void setupRegions( ActionVolume<VolumeAround, TPM>* action,
48 : const Pbc& pbc,
49 : const std::vector<Vector>& positions ) {}
50 : template<typename TPM>
51 : static void parseAtoms( ActionVolume<VolumeAround, TPM>* action,
52 : std::vector<AtomNumber>& atom );
53 : static void calculateNumberInside( const VolumeInput& input,
54 : const VolumeAround& actioninput,
55 : VolumeOutput& output );
56 : #ifdef __PLUMED_HAS_OPENACC
57 : void toACCDevice() const {
58 : #pragma acc enter data copyin(this[0:1],dox,doy,doz,sigma,\
59 : xlow,xhigh,ylow,yhigh,zlow,zhigh,kerneltype)
60 : }
61 : void removeFromACCDevice() const {
62 : #pragma acc exit data delete(kerneltype,zhigh,zlow,yhigh,ylow,xhigh,xlow,\
63 : sigma,doz,doy,dox,this[0:1])
64 : }
65 : #endif //__PLUMED_HAS_OPENACC
66 : };
67 :
68 169 : void VolumeAround::registerKeywords( Keywords& keys ) {
69 338 : keys.setDisplayName("AROUND");
70 169 : keys.add("atoms","ORIGIN","the atom whose vicinity we are interested in examining");
71 169 : keys.addDeprecatedKeyword("ATOM","ORIGIN");
72 169 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
73 169 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
74 169 : keys.add("compulsory","XLOWER","0.0","the lower boundary in x relative to the x coordinate of the atom (0 indicates use full extent of box).");
75 169 : keys.add("compulsory","XUPPER","0.0","the upper boundary in x relative to the x coordinate of the atom (0 indicates use full extent of box).");
76 169 : keys.add("compulsory","YLOWER","0.0","the lower boundary in y relative to the y coordinate of the atom (0 indicates use full extent of box).");
77 169 : keys.add("compulsory","YUPPER","0.0","the upper boundary in y relative to the y coordinate of the atom (0 indicates use full extent of box).");
78 169 : keys.add("compulsory","ZLOWER","0.0","the lower boundary in z relative to the z coordinate of the atom (0 indicates use full extent of box).");
79 169 : keys.add("compulsory","ZUPPER","0.0","the upper boundary in z relative to the z coordinate of the atom (0 indicates use full extent of box).");
80 169 : }
81 :
82 : template<typename TPM>
83 49 : void VolumeAround::parseAtoms( ActionVolume<VolumeAround, TPM>* action, std::vector<AtomNumber>& atom ) {
84 98 : action->parseAtomList("ORIGIN",atom);
85 49 : if( atom.size()==0 ) {
86 0 : action->parseAtomList("ATOM",atom);
87 : }
88 49 : if( atom.size()!=1 ) {
89 0 : action->error("should only be one atom specified");
90 : }
91 49 : action->log.printf(" boundaries for region are calculated based on positions of atom : %d\n",atom[0].serial() );
92 49 : }
93 :
94 : template<typename TPM>
95 49 : void VolumeAround::parseInput( ActionVolume<VolumeAround, TPM>* action ) {
96 98 : action->parse("SIGMA",sigma);
97 : std::string mykerneltype;
98 49 : action->parse("KERNEL",mykerneltype);
99 49 : kerneltype=HistogramBead::getKernelType(mykerneltype);
100 49 : dox=true;
101 49 : action->parse("XLOWER",xlow);
102 49 : action->parse("XUPPER",xhigh);
103 49 : doy=true;
104 49 : action->parse("YLOWER",ylow);
105 49 : action->parse("YUPPER",yhigh);
106 49 : doz=true;
107 49 : action->parse("ZLOWER",zlow);
108 49 : action->parse("ZUPPER",zhigh);
109 49 : if( xlow==0.0 && xhigh==0.0 ) {
110 8 : dox=false;
111 : }
112 49 : if( ylow==0.0 && yhigh==0.0 ) {
113 16 : doy=false;
114 : }
115 49 : if( zlow==0.0 && zhigh==0.0 ) {
116 16 : doz=false;
117 : }
118 49 : if( !dox && !doy && !doz ) {
119 0 : action->error("no subregion defined use XLOWER, XUPPER, YLOWER, YUPPER, ZLOWER, ZUPPER");
120 : }
121 49 : action->log.printf(" boundaries for region (region of interest about atom) : x %f %f, y %f %f, z %f %f \n",xlow,xhigh,ylow,yhigh,zlow,zhigh);
122 49 : }
123 :
124 140716 : void VolumeAround::calculateNumberInside( const VolumeInput& input,
125 : const VolumeAround& actioninput,
126 : VolumeOutput& output ) {
127 : // Setup the histogram bead
128 140716 : HistogramBead bead(actioninput.kerneltype, actioninput.xlow, actioninput.xhigh, actioninput.sigma );
129 :
130 : // Calculate position of atom wrt to origin
131 140716 : Vector fpos=input.pbc.distance( Vector(input.refpos[0][0],input.refpos[0][1],input.refpos[0][2]),
132 140716 : Vector(input.cpos[0],input.cpos[1],input.cpos[2]) );
133 : double xcontr=1.0;
134 140716 : double xder=0.0;
135 140716 : if( actioninput.dox ) {
136 : //bead parameters set in the constructor
137 96716 : xcontr=bead.calculate( fpos[0], xder );
138 : }
139 : double ycontr=1.0;
140 140716 : double yder=0.0;
141 140716 : if( actioninput.doy ) {
142 92288 : bead.set( actioninput.ylow, actioninput.yhigh, actioninput.sigma );
143 92288 : ycontr=bead.calculate( fpos[1], yder );
144 : }
145 : double zcontr=1.0;
146 140716 : double zder=0.0;
147 140716 : if( actioninput.doz ) {
148 92288 : bead.set( actioninput.zlow, actioninput.zhigh, actioninput.sigma );
149 92288 : zcontr=bead.calculate( fpos[2], zder );
150 : }
151 :
152 140716 : output.derivatives[0]=xder*ycontr*zcontr;
153 140716 : output.derivatives[1]=xcontr*yder*zcontr;
154 140716 : output.derivatives[2]=xcontr*ycontr*zder;
155 : // Add derivatives wrt to position of origin atom
156 140716 : output.refders[0][0] = -output.derivatives[0];
157 140716 : output.refders[0][1] = -output.derivatives[1];
158 140716 : output.refders[0][2] = -output.derivatives[2];
159 : // Add virial contribution
160 140716 : output.virial.set( 0, -Tensor(fpos,
161 140716 : Vector(output.derivatives[0], output.derivatives[1], output.derivatives[2])) );
162 140716 : output.values[0] = xcontr*ycontr*zcontr;
163 140716 : }
164 :
165 : }
166 : }
167 : #endif //__PLUMED_volumes_VolumeAround_h
|
