Line data    Source code

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