Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-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 "ActionWithVirtualAtom.h"
      23             : #include "Atoms.h"
      24             : 
      25             : using namespace std;
      26             : 
      27             : namespace PLMD {
      28             : 
      29         150 : void ActionWithVirtualAtom::registerKeywords(Keywords& keys) {
      30         150 :   Action::registerKeywords(keys);
      31         150 :   ActionAtomistic::registerKeywords(keys);
      32         600 :   keys.add("atoms","ATOMS","the list of atoms which are involved the virtual atom's definition");
      33         150 : }
      34             : 
      35         145 : ActionWithVirtualAtom::ActionWithVirtualAtom(const ActionOptions&ao):
      36             :   Action(ao),
      37             :   ActionAtomistic(ao),
      38         290 :   boxDerivatives(3)
      39             : {
      40         145 :   index=atoms.addVirtualAtom(this);
      41         290 :   log.printf("  serial associated to this virtual atom is %u\n",index.serial());
      42         145 : }
      43             : 
      44         290 : ActionWithVirtualAtom::~ActionWithVirtualAtom() {
      45         145 :   atoms.removeVirtualAtom(this);
      46         145 : }
      47             : 
      48        3328 : void ActionWithVirtualAtom::apply() {
      49        3328 :   Vector & f(atoms.forces[index.index()]);
      50       58278 :   for(unsigned i=0; i<getNumberOfAtoms(); i++) modifyForces()[i]=matmul(derivatives[i],f);
      51             :   Tensor & v(modifyVirial());
      52       23296 :   for(unsigned i=0; i<3; i++) v+=boxDerivatives[i]*f[i];
      53        3328 :   f.zero(); // after propagating the force to the atoms used to compute the vatom, we reset this to zero
      54             :   // this is necessary to avoid double counting if then one tries to compute the total force on the c.o.m. of the system.
      55             :   // notice that this is currently done in FIT_TO_TEMPLATE
      56        3328 : }
      57             : 
      58         137 : void ActionWithVirtualAtom::requestAtoms(const std::vector<AtomNumber> & a) {
      59         137 :   ActionAtomistic::requestAtoms(a);
      60         274 :   derivatives.resize(a.size());
      61         137 : }
      62             : 
      63           9 : void ActionWithVirtualAtom::setGradients() {
      64           9 :   gradients.clear();
      65          81 :   for(unsigned i=0; i<getNumberOfAtoms(); i++) {
      66          72 :     AtomNumber an=getAbsoluteIndex(i);
      67             :     // this case if the atom is a virtual one
      68          72 :     if(atoms.isVirtualAtom(an)) {
      69             :       const ActionWithVirtualAtom* a=atoms.getVirtualAtomsAction(an);
      70          36 :       for(const auto & p : a->gradients) {
      71          60 :         gradients[p.first]+=matmul(derivatives[i],p.second);
      72             :       }
      73             :       // this case if the atom is a normal one
      74             :     } else {
      75          60 :       gradients[an]+=derivatives[i];
      76             :     }
      77             :   }
      78           9 : }
      79             : 
      80          36 : void ActionWithVirtualAtom::setBoxDerivatives(const std::vector<Tensor> &d) {
      81          36 :   boxDerivatives=d;
      82          36 :   plumed_assert(d.size()==3);
      83             : // Subtract the trivial part coming from a distorsion applied to the ghost atom first.
      84             : // Notice that this part alone should exactly cancel the already accumulated virial
      85             : // due to forces on this atom.
      86          72 :   Vector pos=atoms.positions[index.index()];
      87         468 :   for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) boxDerivatives[j][i][j]+=pos[i];
      88          36 : }
      89             : 
      90          36 : void ActionWithVirtualAtom::setBoxDerivativesNoPbc() {
      91          36 :   std::vector<Tensor> bd(3);
      92        1116 :   for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) for(unsigned k=0; k<3; k++) {
      93             : // Notice that this expression is very similar to the one used in Colvar::setBoxDerivativesNoPbc().
      94             : // Indeed, we have the negative of a sum over dependent atoms (l) of the external product between positions
      95             : // and derivatives. Notice that this only works only when Pbc have not been used to compute
      96             : // derivatives.
      97        2106 :         for(unsigned l=0; l<getNumberOfAtoms(); l++) {
      98        2268 :           bd[k][i][j]-=getPosition(l)[i]*derivatives[l][j][k];
      99             :         }
     100             :       }
     101          36 :   setBoxDerivatives(bd);
     102          36 : }
     103             : 
     104             : 
     105             : 
     106        3412 : void ActionWithVirtualAtom::setGradientsIfNeeded() {
     107        6824 :   if(isOptionOn("GRADIENTS")) {
     108           9 :     setGradients() ;
     109             :   }
     110        3412 : }
     111             : 
     112        4839 : }