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