LCOV - code coverage report
Current view: top level - generic - Group.cpp (source / functions) Hit Total Coverage
Test: plumed test coverage Lines: 69 71 97.2 %
Date: 2020-11-18 11:20:57 Functions: 9 13 69.2 %

          Line data    Source code
       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-2019 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             : 
      23             : #include "core/ActionRegister.h"
      24             : #include "core/ActionAtomistic.h"
      25             : #include "core/Atoms.h"
      26             : #include "tools/IFile.h"
      27             : #include "tools/Tools.h"
      28             : #include <string>
      29             : #include <vector>
      30             : #include <algorithm>
      31             : 
      32             : using namespace std;
      33             : 
      34             : namespace PLMD {
      35             : namespace generic {
      36             : 
      37             : //+PLUMEDOC GENERIC GROUP
      38             : /*
      39             : Define a group of atoms so that a particular list of atoms can be referenced with a single label
      40             : in definitions of CVs or virtual atoms.
      41             : 
      42             : Atoms can be listed as comma separated numbers (i.e. `1,2,3,10,45,7,9`) , simple positive ranges
      43             : (i.e. `20-40`), ranges with a stride either positive or negative (i.e. `20-40:2` or `80-50:-2`) or as
      44             : comma separated combinations of all the former methods (`1,2,4,5,10-20,21-40:2,80-50:-2`).
      45             : 
      46             : Moreover, lists can be imported from ndx files (GROMACS format). Use `NDX_FILE` to set the name of
      47             : the index file and `NDX_GROUP` to set the name of the group to be imported (default is first one).
      48             : 
      49             : It is also possible to remove atoms from a list and or sort them using keywords `REMOVE`, `SORT`, and `UNIQUE`.
      50             : The flow is the following:
      51             : - If `ATOMS` is present, then take the ordered list of atoms from the `ATOMS` keyword as a starting list.
      52             : - If `NDX_FILE` is present, then append to it the list obtained from the gromacs group.
      53             : - If `REMOVE` is present, then remove the first occurence of each of these atoms from the list.
      54             :   If one tries to remove an atom that was not listed plumed adds a notice in the output.
      55             :   An atom that is present twice in the original list might be removed twice.
      56             : - If `SORT` is present, then the resulting list is sorted by increasing serial number.
      57             : - If `UNIQUE` is present, then the resuling list is sorted by increasing serial number _and_ duplicate elements are removed.
      58             : 
      59             : Notice that this command just creates a shortcut, and does not imply any real calculation.
      60             : So, having a huge group defined does not slow down your calculation in any way.
      61             : It is just convenient to better organize input files. Might be used in combination with
      62             : the \ref INCLUDE command so as to store long group definitions in a separate file.
      63             : 
      64             : 
      65             : \par Examples
      66             : 
      67             : This command create a group of atoms containing atoms 1, 4, 7, 11 and 14 (labeled 'o'), and another containing
      68             : atoms 2, 3, 5, 6, 8, 9, 12, and 13 (labeled 'h'):
      69             : \plumedfile
      70             : o: GROUP ATOMS=1,4,7,11,14
      71             : h: GROUP ATOMS=2,3,5,6,8,9,12,13
      72             : # compute the coordination among the two groups
      73             : c: COORDINATION GROUPA=o GROUPB=h R_0=0.3
      74             : # same could have been obtained without GROUP, just writing:
      75             : # c: COORDINATION GROUPA=1,4,7,11,14 GROUPB=2,3,5,6,8,9,12,13
      76             : 
      77             : # print the coordination on file 'colvar'
      78             : PRINT ARG=c FILE=colvar
      79             : \endplumedfile
      80             : 
      81             : Groups can be conveniently stored in a separate file.
      82             : E.g. one could create a file named `groups.dat` which reads
      83             : \plumedfile
      84             : o: GROUP ATOMS=1,4,7,11,14
      85             : h: GROUP ATOMS=2,3,5,6,8,9,12,13
      86             : \endplumedfile
      87             : and then include it in the main 'plumed.dat' file
      88             : \plumedfile
      89             : INCLUDE FILE=groups.dat
      90             : # compute the coordination among the two groups
      91             : c: COORDINATION GROUPA=o GROUPB=h R_0=0.3
      92             : # print the coordination on file 'colvar'
      93             : PRINT ARG=c FILE=colvar
      94             : \endplumedfile
      95             : The `groups.dat` file could be very long and include lists of thousand atoms without cluttering the main plumed.dat file.
      96             : 
      97             : A GROMACS index file can also be imported
      98             : \plumedfile
      99             : # import group named 'protein' from file index.ndx
     100             : pro: GROUP NDX_FILE=index.ndx NDX_GROUP=protein
     101             : # dump all the atoms of the protein on a trajectory file
     102             : DUMPATOMS ATOMS=pro FILE=traj.gro
     103             : \endplumedfile
     104             : 
     105             : A list can be edited with `REMOVE`. For instance, if you
     106             : are using a water model with three atoms per molecule, you can
     107             : easily construct the list of hydrogens in this manner
     108             : \plumedfile
     109             : # take one atom every three, that is oxygens
     110             : ox: GROUP ATOMS=1-90:3
     111             : # take the remaining atoms, that is hydrogens
     112             : hy: GROUP ATOMS=1-90 REMOVE=ox
     113             : DUMPATOMS ATOMS=ox FILE=ox.gro
     114             : DUMPATOMS ATOMS=hy FILE=hy.gro
     115             : \endplumedfile
     116             : 
     117             : 
     118             : */
     119             : //+ENDPLUMEDOC
     120             : 
     121             : class Group:
     122             :   public ActionAtomistic
     123             : {
     124             : 
     125             : public:
     126             :   explicit Group(const ActionOptions&ao);
     127             :   ~Group();
     128             :   static void registerKeywords( Keywords& keys );
     129           0 :   void calculate() {}
     130           0 :   void apply() {}
     131             : };
     132             : 
     133        6539 : PLUMED_REGISTER_ACTION(Group,"GROUP")
     134             : 
     135          87 : Group::Group(const ActionOptions&ao):
     136             :   Action(ao),
     137          87 :   ActionAtomistic(ao)
     138             : {
     139             :   vector<AtomNumber> atoms;
     140         174 :   parseAtomList("ATOMS",atoms);
     141             :   std::string ndxfile,ndxgroup;
     142         174 :   parse("NDX_FILE",ndxfile);
     143         174 :   parse("NDX_GROUP",ndxgroup);
     144          93 :   if(ndxfile.length()>0 && atoms.size()>0) error("either use explicit atom list or import from index file");
     145         168 :   if(ndxfile.length()==0 && ndxgroup.size()>0) error("NDX_GROUP can be only used is NDX_FILE is also used");
     146             : 
     147          87 :   if(ndxfile.length()>0) {
     148          16 :     if(ndxgroup.size()>0) log<<"  importing group '"+ndxgroup+"'";
     149           1 :     else                  log<<"  importing first group";
     150           6 :     log<<" from index file "<<ndxfile<<"\n";
     151             : 
     152          12 :     IFile ifile;
     153           6 :     ifile.open(ndxfile);
     154             :     std::string line;
     155             :     std::string groupname;
     156             :     bool firstgroup=true;
     157             :     bool groupfound=false;
     158        2468 :     while(ifile.getline(line)) {
     159        2462 :       std::vector<std::string> words=Tools::getWords(line);
     160        2487 :       if(words.size()>=3 && words[0]=="[" && words[2]=="]") {
     161          44 :         if(groupname.length()>0) firstgroup=false;
     162             :         groupname=words[1];
     163          83 :         if(groupname==ndxgroup || ndxgroup.length()==0) groupfound=true;
     164        2505 :       } else if(groupname==ndxgroup || (firstgroup && ndxgroup.length()==0)) {
     165        5928 :         for(unsigned i=0; i<words.size(); i++) {
     166        1888 :           AtomNumber at; Tools::convert(words[i],at);
     167        1888 :           atoms.push_back(at);
     168             :         }
     169             :       }
     170             :     }
     171           6 :     if(!groupfound) error("group has not been found in index file");
     172             :   }
     173             : 
     174             :   std::vector<AtomNumber> remove;
     175         174 :   parseAtomList("REMOVE",remove);
     176          87 :   if(remove.size()>0) {
     177             :     std::vector<AtomNumber> notfound;
     178             :     unsigned k=0;
     179           1 :     log<<"  removing these atoms from the list:";
     180          14 :     for(unsigned i=0; i<remove.size(); i++) {
     181           4 :       const auto it = find(atoms.begin(),atoms.end(),remove[i]);
     182           4 :       if(it!=atoms.end()) {
     183           3 :         if(k%25==0) log<<"\n";
     184           6 :         log<<" "<<(*it).serial();
     185           3 :         k++;
     186             :         atoms.erase(it);
     187           1 :       } else notfound.push_back(remove[i]);
     188             :     }
     189           1 :     log<<"\n";
     190           1 :     if(notfound.size()>0) {
     191           1 :       log<<"  the following atoms were not found:";
     192           5 :       for(unsigned i=0; i<notfound.size(); i++) log<<" "<<notfound[i].serial();
     193           1 :       log<<"\n";
     194             :     }
     195             :   }
     196             : 
     197          87 :   bool sortme=false;
     198         174 :   parseFlag("SORT",sortme);
     199          87 :   if(sortme) {
     200           1 :     log<<"  atoms are sorted\n";
     201             :     sort(atoms.begin(),atoms.end());
     202             :   }
     203          87 :   bool unique=false;
     204         174 :   parseFlag("UNIQUE",unique);
     205          87 :   if(unique) {
     206           1 :     log<<"  sorting atoms and removing duplicates\n";
     207           1 :     Tools::removeDuplicates(atoms);
     208             :   }
     209             : 
     210         174 :   this->atoms.insertGroup(getLabel(),atoms);
     211          87 :   log.printf("  list of atoms:");
     212       52218 :   for(unsigned i=0; i<atoms.size(); i++) {
     213       17348 :     if(i%25==0) log<<"\n";
     214       34696 :     log<<" "<<atoms[i].serial();
     215             :   }
     216          87 :   log.printf("\n");
     217          87 : }
     218             : 
     219          88 : void Group::registerKeywords( Keywords& keys ) {
     220          88 :   Action::registerKeywords( keys );
     221          88 :   ActionAtomistic::registerKeywords( keys );
     222         352 :   keys.add("atoms", "ATOMS", "the numerical indexes for the set of atoms in the group");
     223         352 :   keys.add("atoms", "REMOVE","remove these atoms from the list");
     224         264 :   keys.addFlag("SORT",false,"sort the resulting list");
     225         264 :   keys.addFlag("UNIQUE",false,"sort atoms and remove duplicated ones");
     226         352 :   keys.add("optional", "NDX_FILE", "the name of index file (gromacs syntax)");
     227         352 :   keys.add("optional", "NDX_GROUP", "the name of the group to be imported (gromacs syntax) - first group found is used by default");
     228          88 : }
     229             : 
     230         261 : Group::~Group() {
     231         174 :   atoms.removeGroup(getLabel());
     232         174 : }
     233             : 
     234             : }
     235        4839 : }
     236             : 

Generated by: LCOV version 1.13