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 <array>
      24             : #include <iostream>
      25             : 
      26             : namespace PLMD {
      27             : 
      28       18068 : void ActionWithVirtualAtom::registerKeywords(Keywords& keys) {
      29       18068 :   Action::registerKeywords(keys);
      30       18068 :   ActionAtomistic::registerKeywords(keys);
      31       36136 :   keys.add("atoms","ATOMS","the list of atoms which are involved the virtual atom's definition");
      32       36136 :   keys.addOutputComponent("x","default","the x coordinate of the virtual atom");
      33       36136 :   keys.addOutputComponent("y","default","the y coordinate of the virtual atom");
      34       36136 :   keys.addOutputComponent("z","default","the z coordinate of the virtual atom");
      35       36136 :   keys.addOutputComponent("mass","default","the mass of the virtual atom");
      36       36136 :   keys.addOutputComponent("charge","default","the charge of the virtual atom");
      37       18068 : }
      38             : 
      39       10320 : ActionWithVirtualAtom::ActionWithVirtualAtom(const ActionOptions&ao):
      40             :   Action(ao),
      41             :   ActionAtomistic(ao),
      42       10320 :   ActionWithValue(ao)
      43             : {
      44       30960 :   addComponentWithDerivatives("x"); componentIsNotPeriodic("x");
      45       30960 :   addComponentWithDerivatives("y"); componentIsNotPeriodic("y");
      46       20640 :   addComponentWithDerivatives("z"); componentIsNotPeriodic("z");
      47             :   // Store the derivatives with respect to the virial only even if there are no atoms
      48       41280 :   for(unsigned i=0; i<3; ++i) getPntrToComponent(i)->resizeDerivatives(9);
      49       41280 :   addComponent("mass"); componentIsNotPeriodic("mass"); getPntrToComponent("mass")->isConstant();
      50       30960 :   addComponent("charge"); componentIsNotPeriodic("charge"); getPntrToComponent("charge")->isConstant();
      51       10320 : }
      52             : 
      53       10284 : void ActionWithVirtualAtom::requestAtoms(const std::vector<AtomNumber> & a) {
      54       41136 :   ActionAtomistic::requestAtoms(a); for(unsigned i=0; i<3; ++i) getPntrToComponent(i)->resizeDerivatives(3*a.size()+9);
      55       10284 : }
      56             : 
      57       35011 : void ActionWithVirtualAtom::apply() {
      58       35011 :   if( !checkForForces() ) return ;
      59             : 
      60       22901 :   Value* xval=getPntrToComponent(0);
      61       22901 :   Value* yval=getPntrToComponent(1);
      62       22901 :   Value* zval=getPntrToComponent(2);
      63       22901 :   if( !xval->forcesWereAdded() && !yval->forcesWereAdded() && !zval->forcesWereAdded() ) return ;
      64       22901 :   if( xval->isConstant() && yval->isConstant() && zval->isConstant() ) return;
      65             : 
      66       45643 :   for(unsigned i=0; i<value_depends.size(); ++i) {
      67       22742 :     xpos[value_depends[i]]->hasForce = true;
      68       22742 :     ypos[value_depends[i]]->hasForce = true;
      69       22742 :     zpos[value_depends[i]]->hasForce = true;
      70             :   }
      71             :   unsigned k=0;
      72       22901 :   double xf = xval->inputForce[0];
      73       22901 :   double yf = yval->inputForce[0];
      74       22901 :   double zf = zval->inputForce[0];
      75             : // for(const auto & a : atom_value_ind) {
      76             : //   std::size_t nn = a.first, kk = a.second;
      77             : //   xpos[nn]->inputForce[kk] += xf*xval->data[1+k] + yf*yval->data[1+k] + zf*zval->data[1+k]; k++;
      78             : //   ypos[nn]->inputForce[kk] += xf*xval->data[1+k] + yf*yval->data[1+k] + zf*zval->data[1+k]; k++;
      79             : //   zpos[nn]->inputForce[kk] += xf*xval->data[1+k] + yf*yval->data[1+k] + zf*zval->data[1+k]; k++;
      80             : // }
      81       45643 :   for(const auto & a : atom_value_ind_grouped) {
      82       22742 :     const auto nn=a.first;
      83       22742 :     auto & xp=xpos[nn]->inputForce;
      84       22742 :     auto & yp=ypos[nn]->inputForce;
      85       22742 :     auto & zp=zpos[nn]->inputForce;
      86      333478 :     for(const auto & kk : a.second) {
      87      310736 :       xp[kk] += xf*xval->data[1+k] + yf*yval->data[1+k] + zf*zval->data[1+k]; k++;
      88      310736 :       yp[kk] += xf*xval->data[1+k] + yf*yval->data[1+k] + zf*zval->data[1+k]; k++;
      89      310736 :       zp[kk] += xf*xval->data[1+k] + yf*yval->data[1+k] + zf*zval->data[1+k]; k++;
      90             :     }
      91             :   }
      92             : 
      93             :   std::array<double,9> virial;
      94      229010 :   for(unsigned i=0; i<9; ++i) { virial[i] = xf*xval->data[1+k] + yf*yval->data[1+k] + zf*zval->data[1+k]; k++; }
      95       22901 :   unsigned ind = 0; setForcesOnCell( virial.data(), virial.size(), ind );
      96             :   // The above can be achieved using the two functions below.  The code above that is specialised for the ActionWithVirtualAtom
      97             :   // class runs faster than the general code below.
      98             :   // if( !checkForForces() ) return ;
      99             :   // unsigned mm=0; setForcesOnAtoms( getForcesToApply(), mm );
     100             : }
     101             : 
     102        2143 : void ActionWithVirtualAtom::setBoxDerivatives(const std::vector<Tensor> &d) {
     103        2143 :   plumed_assert(d.size()==3); unsigned nbase = 3*getNumberOfAtoms();
     104        8572 :   for(unsigned i=0; i<3; ++i) {
     105        6429 :     Value* myval = getPntrToComponent(i);
     106       25716 :     for(unsigned j=0; j<3; ++j) {
     107       77148 :       for(unsigned k=0; k<3; ++k) myval->setDerivative( nbase + 3*j + k, d[i][j][k] );
     108             :     }
     109             :   }
     110             : // Subtract the trivial part coming from a distorsion applied to the ghost atom first.
     111             : // Notice that this part alone should exactly cancel the already accumulated virial
     112             : // due to forces on this atom.
     113        8572 :   Vector pos; for(unsigned i=0; i<3; ++i) pos[i] = getPntrToComponent(i)->get();
     114       27859 :   for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) getPntrToComponent(j)->addDerivative( nbase + 3*i + j, pos[i] );
     115        2143 : }
     116             : 
     117        2143 : void ActionWithVirtualAtom::setBoxDerivativesNoPbc() {
     118        2143 :   std::vector<Tensor> bd(3);
     119       85720 :   for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) for(unsigned k=0; k<3; k++) {
     120             : // Notice that this expression is very similar to the one used in Colvar::setBoxDerivativesNoPbc().
     121             : // Indeed, we have the negative of a sum over dependent atoms (l) of the external product between positions
     122             : // and derivatives. Notice that this only works only when Pbc have not been used to compute
     123             : // derivatives.
     124      172314 :         for(unsigned l=0; l<getNumberOfAtoms(); l++) {
     125      114453 :           bd[k][i][j]-=getPosition(l)[i]*getPntrToComponent(k)->getDerivative(3*l+j);
     126             :         }
     127             :       }
     128        2143 :   setBoxDerivatives(bd);
     129        2143 : }
     130             : 
     131             : }