Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-2018 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         122 : void ActionWithVirtualAtom::registerKeywords(Keywords& keys) {
      30         122 :   Action::registerKeywords(keys);
      31         122 :   ActionAtomistic::registerKeywords(keys);
      32         122 :   keys.add("atoms","ATOMS","the list of atoms which are involved the virtual atom's definition");
      33         122 : }
      34             : 
      35         117 : ActionWithVirtualAtom::ActionWithVirtualAtom(const ActionOptions&ao):
      36             :   Action(ao),
      37             :   ActionAtomistic(ao),
      38         117 :   boxDerivatives(3)
      39             : {
      40         117 :   index=atoms.addVirtualAtom(this);
      41         117 :   log.printf("  serial associated to this virtual atom is %u\n",index.serial());
      42         117 : }
      43             : 
      44         234 : ActionWithVirtualAtom::~ActionWithVirtualAtom() {
      45         117 :   atoms.removeVirtualAtom(this);
      46         117 : }
      47             : 
      48        3250 : void ActionWithVirtualAtom::apply() {
      49        3250 :   Vector & f(atoms.forces[index.index()]);
      50        3250 :   for(unsigned i=0; i<getNumberOfAtoms(); i++) modifyForces()[i]=matmul(derivatives[i],f);
      51        3250 :   Tensor & v(modifyVirial());
      52        3250 :   for(unsigned i=0; i<3; i++) v+=boxDerivatives[i]*f[i];
      53        3250 :   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        3250 : }
      57             : 
      58         112 : void ActionWithVirtualAtom::requestAtoms(const std::vector<AtomNumber> & a) {
      59         112 :   ActionAtomistic::requestAtoms(a);
      60         112 :   derivatives.resize(a.size());
      61         112 : }
      62             : 
      63           9 : void ActionWithVirtualAtom::setGradients() {
      64           9 :   gradients.clear();
      65          45 :   for(unsigned i=0; i<getNumberOfAtoms(); i++) {
      66          36 :     AtomNumber an=getAbsoluteIndex(i);
      67             :     // this case if the atom is a virtual one
      68          36 :     if(atoms.isVirtualAtom(an)) {
      69           6 :       const ActionWithVirtualAtom* a=atoms.getVirtualAtomsAction(an);
      70          36 :       for(std::map<AtomNumber,Tensor>::const_iterator p=a->gradients.begin(); p!=a->gradients.end(); ++p) {
      71          30 :         gradients[(*p).first]+=matmul(derivatives[i],(*p).second);
      72             :       }
      73             :       // this case if the atom is a normal one
      74             :     } else {
      75          30 :       gradients[an]+=derivatives[i];
      76             :     }
      77             :   }
      78           9 : }
      79             : 
      80          15 : void ActionWithVirtualAtom::setBoxDerivatives(const std::vector<Tensor> &d) {
      81          15 :   boxDerivatives=d;
      82          15 :   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          15 :   Vector pos=atoms.positions[index.index()];
      87          15 :   for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) boxDerivatives[j][i][j]+=pos[i];
      88          15 : }
      89             : 
      90          15 : void ActionWithVirtualAtom::setBoxDerivativesNoPbc() {
      91          15 :   std::vector<Tensor> bd(3);
      92         420 :   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         972 :         for(unsigned l=0; l<getNumberOfAtoms(); l++) {
      98         567 :           bd[k][i][j]-=getPosition(l)[i]*derivatives[l][j][k];
      99             :         }
     100             :       }
     101          15 :   setBoxDerivatives(bd);
     102          15 : }
     103             : 
     104             : 
     105             : 
     106        3310 : void ActionWithVirtualAtom::setGradientsIfNeeded() {
     107        3310 :   if(isOptionOn("GRADIENTS")) {
     108           9 :     setGradients() ;
     109             :   }
     110        3310 : }
     111             : 
     112        2523 : }