Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-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 "core/ActionRegister.h"
      23             : #include "VectorMultiColvar.h"
      24             : 
      25             : namespace PLMD {
      26             : namespace crystallization {
      27             : 
      28             : //+PLUMEDOC MCOLVAR MOLECULES
      29             : /*
      30             : Calculate the vectors connecting a pair of atoms in order to represent the orientation of a molecule.
      31             : 
      32             : At its simplest this command can be used to calculate the average length of an internal vector in a
      33             : collection of different molecules.  When used in conjunction with MutiColvarFunctions in can be used
      34             : to do a variety of more complex tasks.
      35             : 
      36             : \par Examples
      37             : 
      38             : The following input tells plumed to calculate the distances between two of the atoms in a molecule.
      39             : This is done for the same set of atoms four different molecules and the average separation is then
      40             : calculated.
      41             : 
      42             : \plumedfile
      43             : MOLECULES MOL1=1,2 MOL2=3,4 MOL3=5,6 MOL4=7,8 MEAN LABEL=mm
      44             : PRINT ARG=mm.mean FILE=colvar
      45             : \endplumedfile
      46             : 
      47             : 
      48             : */
      49             : //+ENDPLUMEDOC
      50             : 
      51             : 
      52             : class MoleculeOrientation : public VectorMultiColvar {
      53             : private:
      54             :   unsigned nvectors;
      55             : public:
      56             :   static void registerKeywords( Keywords& keys );
      57             :   explicit MoleculeOrientation( const ActionOptions& ao );
      58             :   AtomNumber getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const override;
      59             :   void calculateVector( multicolvar::AtomValuePack& myatoms ) const override;
      60             :   void normalizeVector( std::vector<double>& vals ) const override;
      61             :   void normalizeVectorDerivatives( MultiValue& myvals ) const override;
      62             : };
      63             : 
      64       10427 : PLUMED_REGISTER_ACTION(MoleculeOrientation,"MOLECULES")
      65             : 
      66           5 : void MoleculeOrientation::registerKeywords( Keywords& keys ) {
      67          10 :   VectorMultiColvar::registerKeywords( keys ); keys.use("MEAN"); keys.use("VMEAN");
      68          10 :   keys.add("numbered","MOL","The numerical indices of the atoms in the molecule. The orientation of the molecule is equal to "
      69             :            "the vector connecting the first two atoms specified.  If a third atom is specified its position "
      70             :            "is used to specify where the molecule is.  If a third atom is not present the molecule is assumed "
      71             :            "to be at the center of the vector connecting the first two atoms.");
      72          10 :   keys.reset_style("MOL","atoms");
      73           5 : }
      74             : 
      75           4 : MoleculeOrientation::MoleculeOrientation( const ActionOptions& ao ):
      76             :   Action(ao),
      77           4 :   VectorMultiColvar(ao)
      78             : {
      79             :   std::vector<AtomNumber> all_atoms;
      80           8 :   readAtomsLikeKeyword("MOL",-1,all_atoms);
      81           4 :   nvectors = std::floor( ablocks.size() / 2 );
      82           4 :   if( ablocks.size()%2!=0 && 2*nvectors+1!=ablocks.size() ) error("number of atoms in molecule specification is wrong.  Should be two or three.");
      83             : 
      84           4 :   if( all_atoms.size()==0 ) error("No atoms were specified");
      85           4 :   setVectorDimensionality( 3*nvectors ); setupMultiColvarBase( all_atoms );
      86             : 
      87           4 :   if( ablocks.size()==2*nvectors+1  ) {
      88           4 :     std::vector<bool> catom_ind(ablocks.size(), false); catom_ind[ablocks.size()-1]=true;
      89           4 :     setAtomsForCentralAtom( catom_ind );
      90             :   }
      91           4 : }
      92             : 
      93           0 : AtomNumber MoleculeOrientation::getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const {
      94             :   plumed_dbg_assert( iatom<atom_lab.size() );
      95           0 :   plumed_assert( atom_lab[iatom].first==0 );
      96           0 :   return ActionAtomistic::getAbsoluteIndex( ablocks[0][atom_lab[iatom].second] );
      97             : }
      98             : 
      99       16584 : void MoleculeOrientation::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
     100       33168 :   for(unsigned i=0; i<nvectors; ++i) {
     101       16584 :     Vector distance; distance=getSeparation( myatoms.getPosition(2*i+0), myatoms.getPosition(2*i+1) );
     102             : 
     103       16584 :     addAtomDerivatives( 2+3*i+0, 2*i+0, Vector(-1.0,0,0), myatoms );
     104       16584 :     addAtomDerivatives( 2+3*i+0, 2*i+1, Vector(+1.0,0,0), myatoms );
     105       16584 :     myatoms.addBoxDerivatives( 2+3*i+0, Tensor(distance,Vector(-1.0,0,0)) );
     106       16584 :     myatoms.addValue( 2+3*i+0, distance[0] );
     107             : 
     108       16584 :     addAtomDerivatives( 2+3*i+1, 2*i+0, Vector(0,-1.0,0), myatoms );
     109       16584 :     addAtomDerivatives( 2+3*i+1, 2*i+1, Vector(0,+1.0,0), myatoms );
     110       16584 :     myatoms.addBoxDerivatives( 2+3*i+1, Tensor(distance,Vector(0,-1.0,0)) );
     111       16584 :     myatoms.addValue( 2+3*i+1, distance[1] );
     112             : 
     113       16584 :     addAtomDerivatives( 2+3*i+2, 2*i+0, Vector(0,0,-1.0), myatoms );
     114       16584 :     addAtomDerivatives( 2+3*i+2, 2*i+1, Vector(0,0,+1.0), myatoms );
     115       16584 :     myatoms.addBoxDerivatives( 2+3*i+2, Tensor(distance,Vector(0,0,-1.0)) );
     116       16584 :     myatoms.addValue( 2+3*i+2, distance[2] );
     117             :   }
     118       16584 : }
     119             : 
     120     1041228 : void MoleculeOrientation::normalizeVector( std::vector<double>& vals ) const {
     121     2082456 :   for(unsigned i=0; i<nvectors; ++i) {
     122             :     double norm=0;
     123     4164912 :     for(unsigned j=0; j<3; ++j) norm += vals[2+3*i+j]*vals[2+3*i+j];
     124     1041228 :     norm = sqrt(norm);
     125             : 
     126     1041228 :     double inorm = 1.0; if( norm>epsilon ) inorm = 1.0 / norm;
     127     4164912 :     for(unsigned j=0; j<3; ++j) vals[2+3*i+j] = inorm*vals[2+3*i+j];
     128             :   }
     129     1041228 : }
     130             : 
     131        4992 : void MoleculeOrientation::normalizeVectorDerivatives( MultiValue& myvals ) const {
     132        4992 :   std::vector<double> weight( nvectors ), wdf( nvectors );
     133        9984 :   for(unsigned ivec=0; ivec<nvectors; ++ivec) {
     134       19968 :     double v=0; for(unsigned jcomp=0; jcomp<3; ++jcomp) v += myvals.get( 2+3*ivec+jcomp )*myvals.get( 2+3*ivec+jcomp );
     135        4992 :     v=sqrt(v); weight[ivec]=1.0; wdf[ivec]=1.0;
     136        4992 :     if( v>epsilon ) { weight[ivec] = 1.0 / v; wdf[ivec] = 1.0 / ( v*v*v ); }
     137             :   }
     138             : 
     139       54672 :   for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
     140       49680 :     unsigned jder=myvals.getActiveIndex(j);
     141       99360 :     for(unsigned ivec=0; ivec<nvectors; ++ivec) {
     142      198720 :       double comp2=0.0; for(unsigned jcomp=0; jcomp<3; ++jcomp) comp2 += myvals.get(2+3*ivec+jcomp)*myvals.getDerivative( 2+3*ivec+jcomp, jder );
     143      198720 :       for(unsigned jcomp=0; jcomp<3; ++jcomp) {
     144      149040 :         myvals.setDerivative( 2+3*ivec+jcomp, jder, weight[ivec]*myvals.getDerivative( 2+3*ivec+jcomp, jder ) - wdf[ivec]*comp2*myvals.get(2+3*ivec+jcomp) );
     145             :       }
     146             :     }
     147             :   }
     148        4992 : }
     149             : 
     150             : }
     151             : }