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