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 "ReferenceAtoms.h"
      23             : #include "core/GenericMolInfo.h"
      24             : #include "tools/OFile.h"
      25             : #include "tools/PDB.h"
      26             : 
      27             : namespace PLMD {
      28             : 
      29         751 : ReferenceAtoms::ReferenceAtoms( const ReferenceConfigurationOptions& ro ):
      30             :   ReferenceConfiguration(ro),
      31         751 :   checks_were_disabled(false) {
      32         751 : }
      33             : 
      34         688 : void ReferenceAtoms::readAtomsFromPDB( const PDB& pdb, const bool allowblocks  ) {
      35         688 :   if( !allowblocks && pdb.getNumberOfAtomBlocks()!=1 ) {
      36           0 :     error("found multi-atom-block pdb format but expecting only one block of atoms");
      37             :   }
      38             : 
      39         688 :   indices.resize(0);
      40         688 :   reference_atoms.resize(0);
      41         688 :   align.resize(0);
      42         688 :   displace.resize(0);
      43         688 :   atom_der_index.resize(0);
      44       10045 :   for(unsigned i=0; i<pdb.size(); ++i) {
      45        9357 :     indices.push_back( pdb.getAtomNumbers()[i] );
      46        9357 :     reference_atoms.push_back( pdb.getPositions()[i] );
      47        9357 :     align.push_back( pdb.getOccupancy()[i] );
      48        9357 :     displace.push_back( pdb.getBeta()[i] );
      49        9357 :     atom_der_index.push_back(i);
      50             :   }
      51         688 : }
      52             : 
      53             : // void ReferenceAtoms::setAtomNumbers( const std::vector<AtomNumber>& numbers ) {
      54             : //   reference_atoms.resize( numbers.size() ); align.resize( numbers.size() );
      55             : //   displace.resize( numbers.size() ); atom_der_index.resize( numbers.size() );
      56             : //   indices.resize( numbers.size() );
      57             : //   for(unsigned i=0; i<numbers.size(); ++i) {
      58             : //     indices[i]=numbers[i]; atom_der_index[i]=i;
      59             : //   }
      60             : // }
      61             : 
      62             : // void ReferenceAtoms::printAtoms( OFile& ofile, const double& lunits ) const {
      63             : //   plumed_assert( indices.size()==reference_atoms.size() && align.size()==reference_atoms.size() && displace.size()==reference_atoms.size() );
      64             : //   for(unsigned i=0; i<reference_atoms.size(); ++i) {
      65             : //     ofile.printf("ATOM  %5d  X   RES  %4u    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
      66             : //                  indices[i].serial(), i,
      67             : //                  lunits*reference_atoms[i][0], lunits*reference_atoms[i][1], lunits*reference_atoms[i][2],
      68             : //                  align[i], displace[i] );
      69             : //   }
      70             : // }
      71             : 
      72             : // bool ReferenceAtoms::parseAtomList( const std::string& key, std::vector<unsigned>& numbers ) {
      73             : //   plumed_assert( numbers.size()==0 );
      74             : //
      75             : //   std::vector<std::string> strings;
      76             : //   if( !parseVector(key,strings,true) ) return false;
      77             : //   Tools::interpretRanges(strings);
      78             : //
      79             : //   numbers.resize( strings.size() );
      80             : //   for(unsigned i=0; i<strings.size(); ++i) {
      81             : //     AtomNumber atom;
      82             : //     if( !Tools::convert(strings[i],atom ) ) error("could not convert " + strings[i] + " into atom number");
      83             : //
      84             : //     bool found=false;
      85             : //     for(unsigned j=0; j<indices.size(); ++j) {
      86             : //       if( atom==indices[j] ) { found=true; numbers[i]=j; break; }
      87             : //     }
      88             : //     if(!found) error("atom labelled " + strings[i] + " is not present in pdb input file");
      89             : //   }
      90             : //   return true;
      91             : // }
      92             : 
      93         455 : void ReferenceAtoms::getAtomRequests( std::vector<AtomNumber>& numbers, bool disable_checks ) {
      94         455 :   singleDomainRequests(numbers,disable_checks);
      95         455 : }
      96             : 
      97         455 : void ReferenceAtoms::singleDomainRequests( std::vector<AtomNumber>& numbers, bool disable_checks ) {
      98         455 :   checks_were_disabled=disable_checks;
      99         455 :   atom_der_index.resize( indices.size() );
     100             : 
     101         455 :   if( numbers.size()==0 ) {
     102        4220 :     for(unsigned i=0; i<indices.size(); ++i) {
     103        4113 :       numbers.push_back( indices[i] );
     104        4113 :       atom_der_index[i]=i;
     105             :     }
     106             :   } else {
     107         348 :     if(!disable_checks) {
     108         348 :       if( numbers.size()!=indices.size() ) {
     109           0 :         error("mismatched numbers of atoms in pdb frames");
     110             :       }
     111             :     }
     112             : 
     113        4812 :     for(unsigned i=0; i<indices.size(); ++i) {
     114             :       bool found=false;
     115        4464 :       if(!disable_checks) {
     116        4464 :         if( indices[i]!=numbers[i] ) {
     117           0 :           error("found mismatched reference atoms in pdb frames");
     118             :         }
     119        4464 :         atom_der_index[i]=i;
     120             :       } else {
     121           0 :         for(unsigned j=0; j<numbers.size(); ++j) {
     122           0 :           if( indices[i]==numbers[j] ) {
     123             :             found=true;
     124           0 :             atom_der_index[i]=j;
     125             :             break;
     126             :           }
     127             :         }
     128             :         if( !found ) {
     129           0 :           atom_der_index[i]=numbers.size();
     130           0 :           numbers.push_back( indices[i] );
     131             :         }
     132             :       }
     133             :     }
     134             :   }
     135         455 : }
     136             : 
     137         494 : void ReferenceAtoms::displaceReferenceAtoms( const double& weight, const std::vector<Vector>& dir ) {
     138             :   plumed_dbg_assert( dir.size()==reference_atoms.size() );
     139         715 :   for(unsigned i=0; i<dir.size(); ++i) {
     140         221 :     reference_atoms[i] += weight*dir.size()*dir[i];
     141             :   }
     142         494 : }
     143             : 
     144             : }