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        7198 : void ActionWithVirtualAtom::registerKeywords(Keywords& keys) {
      28        7198 :   Action::registerKeywords(keys);
      29        7198 :   ActionAtomistic::registerKeywords(keys);
      30       14396 :   keys.add("atoms","ATOMS","the list of atoms which are involved the virtual atom's definition");
      31        7198 : }
      32             : 
      33        7178 : ActionWithVirtualAtom::ActionWithVirtualAtom(const ActionOptions&ao):
      34             :   Action(ao),
      35             :   ActionAtomistic(ao),
      36        7178 :   index(atoms.addVirtualAtom(this)) {
      37        7178 :   log<<"  serial associated to this virtual atom is "<<index.serial()<<"\n";
      38        7178 : }
      39             : 
      40        7178 : ActionWithVirtualAtom::~ActionWithVirtualAtom() {
      41        7178 :   atoms.removeVirtualAtom(this);
      42        7178 : }
      43             : 
      44       14338 : void ActionWithVirtualAtom::apply() {
      45       14338 :   Vector & f(atoms.forces[index.index()]);
      46       94976 :   for(unsigned i=0; i<getNumberOfAtoms(); i++) {
      47       80638 :     modifyForces()[i]=matmul(derivatives[i],f);
      48             :   }
      49             :   Tensor & v(modifyVirial());
      50       57352 :   for(unsigned i=0; i<3; i++) {
      51       43014 :     v+=boxDerivatives[i]*f[i];
      52             :   }
      53       14338 :   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       14338 : }
      57             : 
      58        7167 : void ActionWithVirtualAtom::requestAtoms(const std::vector<AtomNumber> & a) {
      59        7167 :   ActionAtomistic::requestAtoms(a);
      60        7167 :   derivatives.resize(a.size());
      61        7167 : }
      62             : 
      63           9 : void ActionWithVirtualAtom::setGradients() {
      64             :   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             :       const ActionWithVirtualAtom* a=atoms.getVirtualAtomsAction(an);
      70          36 :       for(const auto & p : a->gradients) {
      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         605 : void ActionWithVirtualAtom::setBoxDerivatives(const std::array<Tensor,3> &d) {
      81         605 :   boxDerivatives=d;
      82             : // Subtract the trivial part coming from a distorsion applied to the ghost atom first.
      83             : // Notice that this part alone should exactly cancel the already accumulated virial
      84             : // due to forces on this atom.
      85         605 :   Vector pos=atoms.positions[index.index()];
      86        2420 :   for(unsigned i=0; i<3; i++)
      87        7260 :     for(unsigned j=0; j<3; j++) {
      88        5445 :       boxDerivatives[j][i][j]+=pos[i];
      89             :     }
      90         605 : }
      91             : 
      92         605 : void ActionWithVirtualAtom::setBoxDerivativesNoPbc() {
      93         605 :   std::array<Tensor,3> bd;
      94        2420 :   for(unsigned i=0; i<3; i++)
      95        7260 :     for(unsigned j=0; j<3; j++)
      96       21780 :       for(unsigned k=0; k<3; k++) {
      97             : // Notice that this expression is very similar to the one used in Colvar::setBoxDerivativesNoPbc().
      98             : // Indeed, we have the negative of a sum over dependent atoms (l) of the external product between positions
      99             : // and derivatives. Notice that this only works only when Pbc have not been used to compute
     100             : // derivatives.
     101       16902 :         for(unsigned l=0; l<getNumberOfAtoms(); l++) {
     102         567 :           bd[k][i][j]-=getPosition(l)[i]*derivatives[l][j][k];
     103             :         }
     104             :       }
     105         605 :   setBoxDerivatives(bd);
     106         605 : }
     107             : 
     108             : 
     109             : 
     110       14422 : void ActionWithVirtualAtom::setGradientsIfNeeded() {
     111       28844 :   if(isOptionOn("GRADIENTS")) {
     112           9 :     setGradients() ;
     113             :   }
     114       14422 : }
     115             : 
     116             : }