1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2014-2021 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 "VolumeGradientBase.h"
23 #include "core/PlumedMain.h"
24 #include "core/ActionSet.h"
25 #include "CatomPack.h"
26
27 namespace PLMD {
28 namespace multicolvar {
29
registerKeywords(Keywords & keys)30 void VolumeGradientBase::registerKeywords( Keywords& keys ) {
31 BridgedMultiColvarFunction::registerKeywords( keys );
32 }
33
VolumeGradientBase(const ActionOptions & ao)34 VolumeGradientBase::VolumeGradientBase(const ActionOptions&ao):
35 Action(ao),
36 BridgedMultiColvarFunction(ao)
37 {
38 }
39
requestAtoms(const std::vector<AtomNumber> & atoms)40 void VolumeGradientBase::requestAtoms( const std::vector<AtomNumber>& atoms ) {
41 ActionAtomistic::requestAtoms(atoms); bridgeVariable=3*atoms.size();
42 std::map<std::string,bool> checklabs;
43 for(const auto & p : getDependencies() ) checklabs.insert(std::pair<std::string,bool>(p->getLabel(),false));
44 for(const auto & p : plumed.getActionSet() ) {
45 if( p->getLabel()==getPntrToMultiColvar()->getLabel() ) break;
46 if( checklabs.count(p->getLabel()) ) checklabs[p->getLabel()]=true;
47 }
48 for(const auto & p : checklabs ) {
49 if( !p.second ) error("the input for the virtual atoms used in the input for this action must appear in the input file before the input multicolvar");
50 }
51 addDependency( getPntrToMultiColvar() );
52 tmpforces.resize( 3*atoms.size()+9 );
53 }
54
doJobsRequiredBeforeTaskList()55 void VolumeGradientBase::doJobsRequiredBeforeTaskList() {
56 ActionWithValue::clearDerivatives();
57 retrieveAtoms(); setupRegions();
58 ActionWithVessel::doJobsRequiredBeforeTaskList();
59 }
60
completeTask(const unsigned & curr,MultiValue & invals,MultiValue & outvals) const61 void VolumeGradientBase::completeTask( const unsigned& curr, MultiValue& invals, MultiValue& outvals ) const {
62 if( getPntrToMultiColvar()->isDensity() ) {
63 outvals.setValue(0, 1.0); outvals.setValue(1, 1.0);
64 } else {
65 // Copy derivatives of the colvar and the value of the colvar
66 invals.copyValues( outvals );
67 if( derivativesAreRequired() ) invals.copyDerivatives( outvals );
68 }
69 calculateAllVolumes( curr, outvals );
70 }
71
setNumberInVolume(const unsigned & ivol,const unsigned & curr,const double & weight,const Vector & wdf,const Tensor & virial,const std::vector<Vector> & refders,MultiValue & outvals) const72 void VolumeGradientBase::setNumberInVolume( const unsigned& ivol, const unsigned& curr, const double& weight,
73 const Vector& wdf, const Tensor& virial, const std::vector<Vector>& refders,
74 MultiValue& outvals ) const {
75 MultiColvarBase* mcolv=getPntrToMultiColvar();
76 if( !mcolv->weightHasDerivatives ) {
77 outvals.setValue(ivol, weight );
78 if( derivativesAreRequired() ) {
79 CatomPack catom; mcolv->getCentralAtomPack( 0, curr, catom );
80 for(unsigned i=0; i<catom.getNumberOfAtomsWithDerivatives(); ++i) {
81 unsigned jatom=3*catom.getIndex(i);
82 outvals.addDerivative( ivol, jatom+0, catom.getDerivative(i,0,wdf) );
83 outvals.addDerivative( ivol, jatom+1, catom.getDerivative(i,1,wdf) );
84 outvals.addDerivative( ivol, jatom+2, catom.getDerivative(i,2,wdf) );
85 }
86 unsigned nmder=getPntrToMultiColvar()->getNumberOfDerivatives();
87 for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) outvals.addDerivative( ivol, nmder-9+3*i+j, virial(i,j) );
88 for(unsigned i=0; i<refders.size(); ++i) {
89 unsigned iatom=nmder+3*i;
90
91 outvals.addDerivative( ivol, iatom+0, refders[i][0] );
92 outvals.addDerivative( ivol, iatom+1, refders[i][1] );
93 outvals.addDerivative( ivol, iatom+2, refders[i][2] );
94 }
95 }
96 } else if(ivol==0) {
97 double ww=outvals.get(0); outvals.setValue(ivol,ww*weight);
98 if( derivativesAreRequired() ) {
99 plumed_merror("This needs testing");
100 CatomPack catom; mcolv->getCentralAtomPack( 0, curr, catom );
101 for(unsigned i=0; i<catom.getNumberOfAtomsWithDerivatives(); ++i) {
102 unsigned jatom=3*catom.getIndex(i);
103 outvals.addDerivative( ivol, jatom+0, weight*outvals.getDerivative(ivol,jatom+0) + ww*catom.getDerivative(i,0,wdf) );
104 outvals.addDerivative( ivol, jatom+1, weight*outvals.getDerivative(ivol,jatom+1) + ww*catom.getDerivative(i,1,wdf) );
105 outvals.addDerivative( ivol, jatom+2, weight*outvals.getDerivative(ivol,jatom+2) + ww*catom.getDerivative(i,2,wdf) );
106 }
107 unsigned nmder=getPntrToMultiColvar()->getNumberOfDerivatives();
108 for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) outvals.addDerivative( ivol, nmder-9+3*i+j, ww*virial(i,j) );
109 for(unsigned i=0; i<refders.size(); ++i) {
110 unsigned iatom=nmder+3*i;
111 outvals.addDerivative( ivol, iatom+0, ww*refders[i][0] );
112 outvals.addDerivative( ivol, iatom+1, ww*refders[i][1] );
113 outvals.addDerivative( ivol, iatom+2, ww*refders[i][2] );
114 }
115 }
116 } else {
117 double ww=outvals.get(0); outvals.setValue(ivol,ww*weight);
118 if( derivativesAreRequired() ) {
119 plumed_merror("This needs testing");
120 CatomPack catom; mcolv->getCentralAtomPack( 0, curr, catom );
121 for(unsigned i=0; i<catom.getNumberOfAtomsWithDerivatives(); ++i) {
122 unsigned jatom=3*catom.getIndex(i);
123 outvals.addDerivative( ivol, jatom+0, ww*catom.getDerivative(i,0,wdf) );
124 outvals.addDerivative( ivol, jatom+1, ww*catom.getDerivative(i,1,wdf) );
125 outvals.addDerivative( ivol, jatom+2, ww*catom.getDerivative(i,2,wdf) );
126 }
127 unsigned nmder=getPntrToMultiColvar()->getNumberOfDerivatives();
128 for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) outvals.addDerivative( ivol, nmder-9+3*i+j, ww*virial(i,j) );
129 for(unsigned i=0; i<refders.size(); ++i) {
130 unsigned iatom=nmder+3*i;
131 outvals.addDerivative( ivol, iatom+0, ww*refders[i][0] );
132 outvals.addDerivative( ivol, iatom+1, ww*refders[i][1] );
133 outvals.addDerivative( ivol, iatom+2, ww*refders[i][2] );
134 }
135 }
136 }
137 }
138
addBridgeForces(const std::vector<double> & bb)139 void VolumeGradientBase::addBridgeForces( const std::vector<double>& bb ) {
140 plumed_dbg_assert( bb.size()==tmpforces.size()-9 );
141 // Forces on local atoms
142 for(unsigned i=0; i<bb.size(); ++i) tmpforces[i]=bb[i];
143 // Virial contribution is zero
144 for(unsigned i=bb.size(); i<bb.size()+9; ++i) tmpforces[i]=0.0;
145 setForcesOnAtoms( tmpforces, 0 );
146 }
147
148 }
149 }
150