Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2017 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 "ActionVolume.h"
23 : #include "core/ActionToPutData.h"
24 : #include "core/PlumedMain.h"
25 :
26 : namespace PLMD {
27 : namespace volumes {
28 :
29 217 : void ActionVolume::registerKeywords( Keywords& keys ) {
30 217 : ActionWithVector::registerKeywords( keys );
31 217 : keys.add("atoms","ATOMS","the group of atoms that you would like to investigate");
32 217 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
33 217 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
34 217 : keys.addFlag("OUTSIDE",false,"calculate quantities for colvars that are on atoms outside the region of interest");
35 434 : keys.setValueDescription("scalar/vector","vector of numbers between 0 and 1 that measure the degree to which each atom is within the volume of interest");
36 217 : }
37 :
38 59 : ActionVolume::ActionVolume(const ActionOptions&ao):
39 : Action(ao),
40 59 : ActionWithVector(ao) {
41 : std::vector<AtomNumber> atoms;
42 118 : parseAtomList("ATOMS",atoms);
43 59 : if( atoms.size()==0 ) {
44 0 : error("no atoms were specified");
45 : }
46 59 : log.printf(" examining positions of atoms ");
47 33354 : for(unsigned i=0; i<atoms.size(); ++i) {
48 33295 : log.printf(" %d", atoms[i].serial() );
49 : }
50 59 : log.printf("\n");
51 59 : ActionAtomistic::requestAtoms( atoms );
52 :
53 59 : parseFlag("OUTSIDE",not_in);
54 59 : sigma=0.0;
55 59 : if( keywords.exists("SIGMA") ) {
56 96 : parse("SIGMA",sigma);
57 : }
58 59 : if( keywords.exists("KERNEL") ) {
59 118 : parse("KERNEL",kerneltype);
60 : }
61 :
62 59 : if( atoms.size()==1 ) {
63 2 : ActionWithValue::addValueWithDerivatives();
64 : } else {
65 58 : std::vector<unsigned> shape(1);
66 58 : shape[0]=atoms.size();
67 58 : ActionWithValue::addValue( shape );
68 : }
69 59 : setNotPeriodic();
70 59 : getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
71 59 : }
72 :
73 104 : bool ActionVolume::isInSubChain( unsigned& nder ) {
74 104 : nder = 0;
75 104 : getFirstActionInChain()->getNumberOfStreamedDerivatives( nder, getPntrToComponent(0) );
76 104 : nder = nder - getNumberOfDerivatives();
77 104 : return true;
78 : }
79 :
80 59 : void ActionVolume::requestAtoms( const std::vector<AtomNumber> & a ) {
81 59 : std::vector<AtomNumber> all_atoms( getAbsoluteIndexes() );
82 128 : for(unsigned i=0; i<a.size(); ++i) {
83 69 : all_atoms.push_back( a[i] );
84 : }
85 59 : ActionAtomistic::requestAtoms( all_atoms );
86 59 : if( getPntrToComponent(0)->getRank()==0 ) {
87 1 : getPntrToComponent(0)->resizeDerivatives( 3*getNumberOfAtoms()+9 );
88 : }
89 59 : }
90 :
91 863 : void ActionVolume::areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions ) {
92 863 : task_reducing_actions.push_back(this);
93 863 : }
94 :
95 533 : void ActionVolume::getNumberOfTasks( unsigned& ntasks ) {
96 533 : setupRegions();
97 533 : ActionWithVector::getNumberOfTasks( ntasks );
98 533 : }
99 :
100 464761 : int ActionVolume::checkTaskStatus( const unsigned& taskno, int& flag ) const {
101 464761 : unsigned nref=getNumberOfAtoms()-getConstPntrToComponent(0)->getShape()[0];
102 464761 : Vector wdf;
103 464761 : Tensor vir;
104 464761 : std::vector<Vector> refders( nref );
105 464761 : double weight=calculateNumberInside( ActionAtomistic::getPosition(taskno), wdf, vir, refders );
106 464761 : if( not_in ) {
107 457 : weight = 1.0 - weight;
108 : }
109 464761 : if( weight>epsilon ) {
110 26062 : return 1;
111 : }
112 : return 0;
113 : }
114 :
115 2081 : void ActionVolume::calculate() {
116 2081 : if( actionInChain() ) {
117 : return;
118 : }
119 1750 : if( getPntrToComponent(0)->getRank()==0 ) {
120 1560 : setupRegions();
121 1560 : unsigned nref = getNumberOfAtoms() - 1;
122 1560 : Vector wdf;
123 1560 : Tensor vir;
124 1560 : std::vector<Vector> refders( nref );
125 1560 : double weight=calculateNumberInside( ActionAtomistic::getPosition(0), wdf, vir, refders );
126 1560 : if( not_in ) {
127 0 : weight = 1.0 - weight;
128 0 : wdf *= -1.;
129 0 : vir *=-1;
130 0 : for(unsigned i=0; i<refders.size(); ++i) {
131 0 : refders[i]*=-1;
132 : }
133 : }
134 : // Atom position
135 1560 : Value* v = getPntrToComponent(0);
136 : v->set( weight );
137 6240 : for(unsigned i=0; i<3; ++i ) {
138 4680 : v->addDerivative( i, wdf[i] );
139 : }
140 : // Add derivatives with respect to reference positions
141 7800 : for(unsigned i=0; i<refders.size(); ++i) {
142 24960 : for(unsigned j=0; j<3; ++j ) {
143 18720 : v->addDerivative( 3 + 3*i + j, refders[i][j] );
144 : }
145 : }
146 : // Add virial
147 : unsigned vbase = 3*getNumberOfAtoms();
148 6240 : for(unsigned i=0; i<3; ++i)
149 18720 : for(unsigned j=0; j<3; ++j) {
150 14040 : v->addDerivative( vbase + 3*i + j, vir(i,j) );
151 : }
152 : } else {
153 190 : runAllTasks();
154 : }
155 : }
156 :
157 62937 : void ActionVolume::performTask( const unsigned& curr, MultiValue& outvals ) const {
158 62937 : unsigned nref=getNumberOfAtoms()-getConstPntrToComponent(0)->getShape()[0];
159 62937 : Vector wdf;
160 62937 : Tensor vir;
161 62937 : std::vector<Vector> refders( nref );
162 62937 : double weight=calculateNumberInside( ActionAtomistic::getPosition(curr), wdf, vir, refders );
163 :
164 62937 : if( not_in ) {
165 4000 : weight = 1.0 - weight;
166 4000 : wdf *= -1.;
167 4000 : vir *=-1;
168 8000 : for(unsigned i=0; i<refders.size(); ++i) {
169 4000 : refders[i]*=-1;
170 : }
171 : }
172 62937 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
173 62937 : outvals.setValue( ostrn, weight );
174 :
175 62937 : if( doNotCalculateDerivatives() ) {
176 : return;
177 : }
178 :
179 : // Atom position
180 146144 : for(unsigned i=0; i<3; ++i ) {
181 109608 : outvals.addDerivative( ostrn, 3*curr+i, wdf[i] );
182 109608 : outvals.updateIndex( ostrn, 3*curr+i );
183 : }
184 : // Add derivatives with respect to reference positions
185 36536 : unsigned vbase = 3*(getNumberOfAtoms()-nref);
186 76572 : for(unsigned i=0; i<refders.size(); ++i) {
187 160144 : for(unsigned j=0; j<3; ++j ) {
188 120108 : outvals.addDerivative( ostrn, vbase, refders[i][j] );
189 120108 : outvals.updateIndex( ostrn, vbase );
190 120108 : vbase++;
191 : }
192 : }
193 : // Add virial
194 146144 : for(unsigned i=0; i<3; ++i) {
195 438432 : for(unsigned j=0; j<3; ++j) {
196 328824 : outvals.addDerivative( ostrn, vbase, vir(i,j) );
197 328824 : outvals.updateIndex( ostrn, vbase );
198 328824 : vbase++;
199 : }
200 : }
201 : }
202 :
203 : }
204 : }
|
