Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2019 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 "core/ActionRegister.h"
23 : #include "GridPrintingBase.h"
24 : #include "tools/OFile.h"
25 :
26 : //+PLUMEDOC GRIDANALYSIS DUMPCUBE
27 : /*
28 : Output a three dimensional grid using the Gaussian cube file format.
29 :
30 : Suppose you have calculated the value of a function on a three dimensional grid.
31 : This function might be a \ref HISTOGRAM or it might be a free energy energy surface
32 : that was calculated from this histogram by using \ref CONVERT_TO_FES. Alternatively,
33 : your function might be a phase-field that was calculated using \ref MULTICOLVARDENS.
34 : Whatever the function is, however, you obviously cannot show it using a typical contour
35 : plotting program such as gnuplot as you have three input variables.
36 :
37 : Tools like VMD have nice features for plotting these types of three dimensional functions
38 : but typically you are required to use a Gaussian cube file format to input the data. This
39 : action thus allows you to output a function evaluated on a grid to a Gaussian cube file format.
40 :
41 : \par Examples
42 :
43 : The input below can be used to postprocess a trajectory. A histogram as a function of the distance
44 : between atoms 1 and 2, the distance between atom 1 and 3 and the angle between the vector connecting atoms 1 and
45 : 2 and 1 and 3 is computed using kernel density estimation. Once all the data contained in the trajectory has been read in and
46 : all the kernels have been added the resulting histogram is output to a file called histoA1.cube. This file has the
47 : Gaussian cube file format. The histogram can thus be visualized using tools such as VMD.
48 :
49 : \plumedfile
50 : x1: DISTANCE ATOMS=1,2
51 : x2: DISTANCE ATOMS=1,3
52 : x3: ANGLE ATOMS=1,2,3
53 :
54 : hA1: HISTOGRAM ARG=x1,x2,x3 GRID_MIN=0.0,0.0,0.0 GRID_MAX=3.0,3.0,3.0 GRID_BIN=10,10,10 BANDWIDTH=1.0,1.0,1.0
55 : DUMPCUBE GRID=hA1 FILE=histoA1.cube
56 : \endplumedfile
57 :
58 : */
59 : //+ENDPLUMEDOC
60 :
61 : namespace PLMD {
62 : namespace gridtools {
63 :
64 12 : class DumpCube : public GridPrintingBase {
65 : private:
66 : unsigned mycomp;
67 : public:
68 : static void registerKeywords( Keywords& keys );
69 : explicit DumpCube(const ActionOptions&ao);
70 : void printGrid( OFile& ofile ) const ;
71 : };
72 :
73 4827 : PLUMED_REGISTER_ACTION(DumpCube,"DUMPCUBE")
74 :
75 7 : void DumpCube::registerKeywords( Keywords& keys ) {
76 7 : GridPrintingBase::registerKeywords( keys );
77 7 : keys.add("optional","COMPONENT","if your input is a vector field use this to specifiy the component of the input vector field for which you wish to output");
78 7 : }
79 :
80 6 : DumpCube::DumpCube(const ActionOptions&ao):
81 : Action(ao),
82 6 : GridPrintingBase(ao)
83 : {
84 6 : fmt = fmt + " ";
85 6 : if( ingrid->getType()!="flat" ) error("cannot dump grid of type " + ingrid->getType() + " using DUMPCUBE");
86 6 : if( ingrid->getDimension()!=3 ) error("cannot print cube file if grid does not contain three dimensional data");
87 :
88 6 : if( ingrid->getNumberOfComponents()==1 ) {
89 6 : mycomp=0;
90 : } else {
91 0 : int tcomp=-1; parse("COMPONENT",tcomp);
92 0 : if( tcomp<0 ) error("component of vector field was not specified - use COMPONENT keyword");
93 0 : mycomp=tcomp*(1+ingrid->getDimension()); if( ingrid->noDerivatives() ) mycomp=tcomp;
94 0 : log.printf(" using %dth component of grid \n",tcomp );
95 : }
96 :
97 6 : checkRead();
98 6 : }
99 :
100 8 : void DumpCube::printGrid( OFile& ofile ) const {
101 8 : double lunit = ingrid->getCubeUnits();
102 :
103 8 : ofile.printf("PLUMED CUBE FILE\n");
104 8 : ofile.printf("OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n");
105 : // Number of atoms followed by position of origin (origin set so that center of grid is in center of cell)
106 8 : std::string ostr = "%d " + fmt + fmt + fmt + "\n";
107 8 : ofile.printf(ostr.c_str(),1,-0.5*lunit*ingrid->getGridExtent(0),-0.5*lunit*ingrid->getGridExtent(1),-0.5*lunit*ingrid->getGridExtent(2));
108 8 : ofile.printf(ostr.c_str(),ingrid->getNbin()[0],lunit*ingrid->getGridSpacing()[0],0.0,0.0); // Number of bins in each direction followed by
109 8 : ofile.printf(ostr.c_str(),ingrid->getNbin()[1],0.0,lunit*ingrid->getGridSpacing()[1],0.0); // shape of voxel
110 8 : ofile.printf(ostr.c_str(),ingrid->getNbin()[2],0.0,0.0,lunit*ingrid->getGridSpacing()[2]);
111 8 : ofile.printf(ostr.c_str(),1,0.0,0.0,0.0); // Fake atom otherwise VMD doesn't work
112 16 : std::vector<unsigned> pp(3); std::vector<unsigned> nbin( ingrid->getNbin() );
113 112 : for(pp[0]=0; pp[0]<nbin[0]; ++pp[0]) {
114 1800 : for(pp[1]=0; pp[1]<nbin[1]; ++pp[1]) {
115 40184 : for(pp[2]=0; pp[2]<nbin[2]; ++pp[2]) {
116 38488 : ofile.printf(fmt.c_str(), ingrid->getGridElement( ingrid->getIndex(pp), mycomp ) );
117 38488 : if(pp[2]%6==5) ofile.printf("\n");
118 : }
119 1696 : ofile.printf("\n");
120 : }
121 8 : }
122 8 : }
123 :
124 : }
125 4821 : }
|
