Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-2017 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 "ActionVolume.h"
      23             : #include "core/ActionToPutData.h"
      24             : #include "core/PlumedMain.h"
      25             : 
      26             : namespace PLMD {
      27             : namespace volumes {
      28             : 
      29         217 : void ActionVolume::registerKeywords( Keywords& keys ) {
      30         217 :   ActionWithVector::registerKeywords( keys );
      31         434 :   keys.add("atoms","ATOMS","the group of atoms that you would like to investigate");
      32         434 :   keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
      33         434 :   keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
      34         434 :   keys.addFlag("OUTSIDE",false,"calculate quantities for colvars that are on atoms outside the region of interest");
      35         217 :   keys.setValueDescription("vector of numbers between 0 and 1 that measure the degree to which each atom is within the volume of interest");
      36         217 : }
      37             : 
      38          59 : ActionVolume::ActionVolume(const ActionOptions&ao):
      39             :   Action(ao),
      40          59 :   ActionWithVector(ao)
      41             : {
      42         118 :   std::vector<AtomNumber> atoms; parseAtomList("ATOMS",atoms);
      43          59 :   if( atoms.size()==0 ) error("no atoms were specified");
      44          59 :   log.printf("  examining positions of atoms ");
      45       33354 :   for(unsigned i=0; i<atoms.size(); ++i) log.printf(" %d", atoms[i].serial() );
      46          59 :   log.printf("\n"); ActionAtomistic::requestAtoms( atoms );
      47             : 
      48          59 :   parseFlag("OUTSIDE",not_in); sigma=0.0;
      49         166 :   if( keywords.exists("SIGMA") ) parse("SIGMA",sigma);
      50         177 :   if( keywords.exists("KERNEL") ) parse("KERNEL",kerneltype);
      51             : 
      52          60 :   if( atoms.size()==1 ) ActionWithValue::addValueWithDerivatives();
      53          58 :   else { std::vector<unsigned> shape(1); shape[0]=atoms.size(); ActionWithValue::addValue( shape ); }
      54          59 :   setNotPeriodic(); getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
      55          59 : }
      56             : 
      57         104 : bool ActionVolume::isInSubChain( unsigned& nder ) {
      58         104 :   nder = 0; getFirstActionInChain()->getNumberOfStreamedDerivatives( nder, getPntrToComponent(0) );
      59         104 :   nder = nder - getNumberOfDerivatives();
      60         104 :   return true;
      61             : }
      62             : 
      63          59 : void ActionVolume::requestAtoms( const std::vector<AtomNumber> & a ) {
      64          59 :   std::vector<AtomNumber> all_atoms( getAbsoluteIndexes() );
      65         128 :   for(unsigned i=0; i<a.size(); ++i) all_atoms.push_back( a[i] );
      66          59 :   ActionAtomistic::requestAtoms( all_atoms );
      67          59 :   if( getPntrToComponent(0)->getRank()==0 ) getPntrToComponent(0)->resizeDerivatives( 3*getNumberOfAtoms()+9 );
      68          59 : }
      69             : 
      70         863 : void ActionVolume::areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions ) {
      71         863 :   task_reducing_actions.push_back(this);
      72         863 : }
      73             : 
      74         533 : void ActionVolume::getNumberOfTasks( unsigned& ntasks ) {
      75         533 :   setupRegions(); ActionWithVector::getNumberOfTasks( ntasks );
      76         533 : }
      77             : 
      78      464761 : int ActionVolume::checkTaskStatus( const unsigned& taskno, int& flag ) const {
      79      464761 :   unsigned nref=getNumberOfAtoms()-getConstPntrToComponent(0)->getShape()[0];
      80      464761 :   Vector wdf; Tensor vir; std::vector<Vector> refders( nref );
      81      464761 :   double weight=calculateNumberInside( ActionAtomistic::getPosition(taskno), wdf, vir, refders );
      82      464761 :   if( not_in ) weight = 1.0 - weight;
      83      464761 :   if( weight>epsilon ) return 1;
      84             :   return 0;
      85             : }
      86             : 
      87        2081 : void ActionVolume::calculate() {
      88        2081 :   if( actionInChain() ) return;
      89        1750 :   if( getPntrToComponent(0)->getRank()==0 ) {
      90        1560 :     setupRegions(); unsigned nref = getNumberOfAtoms() - 1;
      91        1560 :     Vector wdf; Tensor vir; std::vector<Vector> refders( nref );
      92        1560 :     double weight=calculateNumberInside( ActionAtomistic::getPosition(0), wdf, vir, refders );
      93        1560 :     if( not_in ) {
      94           0 :       weight = 1.0 - weight; wdf *= -1.; vir *=-1;
      95           0 :       for(unsigned i=0; i<refders.size(); ++i) refders[i]*=-1;
      96             :     }
      97             :     // Atom position
      98        1560 :     Value* v = getPntrToComponent(0); v->set( weight );
      99        6240 :     for(unsigned i=0; i<3; ++i ) v->addDerivative( i, wdf[i] );
     100             :     // Add derivatives with respect to reference positions
     101        7800 :     for(unsigned i=0; i<refders.size(); ++i) {
     102       24960 :       for(unsigned j=0; j<3; ++j ) v->addDerivative( 3 + 3*i + j, refders[i][j] );
     103             :     }
     104             :     // Add virial
     105             :     unsigned vbase = 3*getNumberOfAtoms();
     106       20280 :     for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) v->addDerivative( vbase + 3*i + j, vir(i,j) );
     107         190 :   } else runAllTasks();
     108             : }
     109             : 
     110       62937 : void ActionVolume::performTask( const unsigned& curr, MultiValue& outvals ) const {
     111       62937 :   unsigned nref=getNumberOfAtoms()-getConstPntrToComponent(0)->getShape()[0];
     112       62937 :   Vector wdf; Tensor vir; std::vector<Vector> refders( nref );
     113       62937 :   double weight=calculateNumberInside( ActionAtomistic::getPosition(curr), wdf, vir, refders );
     114             : 
     115       62937 :   if( not_in ) {
     116        4000 :     weight = 1.0 - weight; wdf *= -1.; vir *=-1;
     117        8000 :     for(unsigned i=0; i<refders.size(); ++i) refders[i]*=-1;
     118             :   }
     119       62937 :   unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
     120       62937 :   outvals.setValue( ostrn, weight );
     121             : 
     122       62937 :   if( doNotCalculateDerivatives() ) return;
     123             : 
     124             :   // Atom position
     125      146144 :   for(unsigned i=0; i<3; ++i ) { outvals.addDerivative( ostrn, 3*curr+i, wdf[i] ); outvals.updateIndex( ostrn, 3*curr+i ); }
     126             :   // Add derivatives with respect to reference positions
     127       36536 :   unsigned vbase = 3*(getNumberOfAtoms()-nref);
     128       76572 :   for(unsigned i=0; i<refders.size(); ++i) {
     129      160144 :     for(unsigned j=0; j<3; ++j ) { outvals.addDerivative( ostrn, vbase, refders[i][j] ); outvals.updateIndex( ostrn, vbase ); vbase++; }
     130             :   }
     131             :   // Add virial
     132      146144 :   for(unsigned i=0; i<3; ++i) {
     133      438432 :     for(unsigned j=0; j<3; ++j) { outvals.addDerivative( ostrn, vbase, vir(i,j) ); outvals.updateIndex( ostrn, vbase ); vbase++; }
     134             :   }
     135             : }
     136             : 
     137             : }
     138             : }